diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index f7260e7..43d2a7b 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -700,7 +700,7 @@ class Spectra(JROData): for pair in pairsList: if pair not in self.pairsList: raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair) - pairsIndexList.append(self.pairsList.index(pair)) + pairsIndexList.append(self.pairsList.index(pair)) for i in range(len(pairsIndexList)): pair = self.pairsList[pairsIndexList[i]] ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0) diff --git a/schainpy/model/graphics/jroplot_data.py b/schainpy/model/graphics/jroplot_data.py index 6ad65ac..8316adf 100644 --- a/schainpy/model/graphics/jroplot_data.py +++ b/schainpy/model/graphics/jroplot_data.py @@ -1,32 +1,33 @@ import os -import zmq import time -import numpy +import glob import datetime -import numpy as np +from multiprocessing import Process + +import zmq +import numpy import matplotlib -import glob -matplotlib.use('TkAgg') import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import make_axes_locatable -from matplotlib.ticker import FuncFormatter, LinearLocator -from multiprocessing import Process +from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator from schainpy.model.proc.jroproc_base import Operation - -plt.ion() +from schainpy.utils import log func = lambda x, pos: ('%s') %(datetime.datetime.fromtimestamp(x).strftime('%H:%M')) -fromtimestamp = lambda x, mintime : (datetime.datetime.utcfromtimestamp(mintime).replace(hour=(x + 5), minute=0) - d1970).total_seconds() +d1970 = datetime.datetime(1970, 1, 1) -d1970 = datetime.datetime(1970,1,1) class PlotData(Operation, Process): + ''' + Base class for Schain plotting operations + ''' CODE = 'Figure' colormap = 'jro' + bgcolor = 'white' CONFLATE = False __MAXNUMX = 80 __missing = 1E30 @@ -37,54 +38,143 @@ class PlotData(Operation, Process): Process.__init__(self) self.kwargs['code'] = self.CODE self.mp = False - self.dataOut = None - self.isConfig = False - self.figure = None + self.data = None + self.isConfig = False + self.figures = [] self.axes = [] + self.cb_axes = [] self.localtime = kwargs.pop('localtime', True) self.show = kwargs.get('show', True) self.save = kwargs.get('save', False) self.colormap = kwargs.get('colormap', self.colormap) self.colormap_coh = kwargs.get('colormap_coh', 'jet') self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r') - self.showprofile = kwargs.get('showprofile', True) - self.title = kwargs.get('wintitle', '') + self.colormaps = kwargs.get('colormaps', None) + self.bgcolor = kwargs.get('bgcolor', self.bgcolor) + self.showprofile = kwargs.get('showprofile', False) + self.title = kwargs.get('wintitle', self.CODE.upper()) + self.cb_label = kwargs.get('cb_label', None) + self.cb_labels = kwargs.get('cb_labels', None) self.xaxis = kwargs.get('xaxis', 'frequency') self.zmin = kwargs.get('zmin', None) self.zmax = kwargs.get('zmax', None) + self.zlimits = kwargs.get('zlimits', None) self.xmin = kwargs.get('xmin', None) + if self.xmin is not None: + self.xmin += 5 self.xmax = kwargs.get('xmax', None) self.xrange = kwargs.get('xrange', 24) self.ymin = kwargs.get('ymin', None) self.ymax = kwargs.get('ymax', None) - self.__MAXNUMY = kwargs.get('decimation', 5000) - self.throttle_value = 5 - self.times = [] - #self.interactive = self.kwargs['parent'] + self.xlabel = kwargs.get('xlabel', None) + self.__MAXNUMY = kwargs.get('decimation', 100) + self.showSNR = kwargs.get('showSNR', False) + self.oneFigure = kwargs.get('oneFigure', True) + self.width = kwargs.get('width', None) + self.height = kwargs.get('height', None) + self.colorbar = kwargs.get('colorbar', True) + self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1]) + self.titles = ['' for __ in range(16)] + + def __setup(self): + ''' + Common setup for all figures, here figures and axes are created + ''' + + self.setup() + + if self.width is None: + self.width = 8 + self.figures = [] + self.axes = [] + self.cb_axes = [] + self.pf_axes = [] + self.cmaps = [] + + size = '15%' if self.ncols==1 else '30%' + pad = '4%' if self.ncols==1 else '8%' + + if self.oneFigure: + if self.height is None: + self.height = 1.4*self.nrows + 1 + fig = plt.figure(figsize=(self.width, self.height), + edgecolor='k', + facecolor='w') + self.figures.append(fig) + for n in range(self.nplots): + ax = fig.add_subplot(self.nrows, self.ncols, n+1) + ax.tick_params(labelsize=8) + ax.firsttime = True + self.axes.append(ax) + if self.showprofile: + cax = self.__add_axes(ax, size=size, pad=pad) + cax.tick_params(labelsize=8) + self.pf_axes.append(cax) + else: + if self.height is None: + self.height = 3 + for n in range(self.nplots): + fig = plt.figure(figsize=(self.width, self.height), + edgecolor='k', + facecolor='w') + ax = fig.add_subplot(1, 1, 1) + ax.tick_params(labelsize=8) + ax.firsttime = True + self.figures.append(fig) + self.axes.append(ax) + if self.showprofile: + cax = self.__add_axes(ax, size=size, pad=pad) + cax.tick_params(labelsize=8) + self.pf_axes.append(cax) + + for n in range(self.nrows): + if self.colormaps is not None: + cmap = plt.get_cmap(self.colormaps[n]) + else: + cmap = plt.get_cmap(self.colormap) + cmap.set_bad(self.bgcolor, 1.) + self.cmaps.append(cmap) + + def __add_axes(self, ax, size='30%', pad='8%'): ''' - this new parameter is created to plot data from varius channels at different figures - 1. crear una lista de figuras donde se puedan plotear las figuras, - 2. dar las opciones de configuracion a cada figura, estas opciones son iguales para ambas figuras - 3. probar? + Add new axes to the given figure ''' - self.ind_plt_ch = kwargs.get('ind_plt_ch', False) - self.figurelist = None + divider = make_axes_locatable(ax) + nax = divider.new_horizontal(size=size, pad=pad) + ax.figure.add_axes(nax) + return nax - def fill_gaps(self, x_buffer, y_buffer, z_buffer): + def setup(self): + ''' + This method should be implemented in the child class, the following + attributes should be set: + + self.nrows: number of rows + self.ncols: number of cols + self.nplots: number of plots (channels or pairs) + self.ylabel: label for Y axes + self.titles: list of axes title + + ''' + raise(NotImplementedError, 'Implement this method in child class') + def fill_gaps(self, x_buffer, y_buffer, z_buffer): + ''' + Create a masked array for missing data + ''' if x_buffer.shape[0] < 2: return x_buffer, y_buffer, z_buffer deltas = x_buffer[1:] - x_buffer[0:-1] - x_median = np.median(deltas) + x_median = numpy.median(deltas) - index = np.where(deltas > 5*x_median) + index = numpy.where(deltas > 5*x_median) if len(index[0]) != 0: z_buffer[::, index[0], ::] = self.__missing - z_buffer = np.ma.masked_inside(z_buffer, + z_buffer = numpy.ma.masked_inside(z_buffer, 0.99*self.__missing, 1.01*self.__missing) @@ -99,110 +189,117 @@ class PlotData(Operation, Process): x = self.x y = self.y[::dy] z = self.z[::, ::, ::dy] - + return x, y, z - ''' - JM: - elimana las otras imagenes generadas debido a que lso workers no llegan en orden y le pueden - poner otro tiempo a la figura q no necesariamente es el ultimo. - Solo se realiza cuando termina la imagen. - Problemas: + def format(self): + ''' + Set min and max values, labels, ticks and titles + ''' - File "/home/ci-81/workspace/schainv2.3/schainpy/model/graphics/jroplot_data.py", line 145, in __plot - for n, eachfigure in enumerate(self.figurelist): - TypeError: 'NoneType' object is not iterable + if self.xmin is None: + xmin = self.min_time + else: + if self.xaxis is 'time': + dt = datetime.datetime.fromtimestamp(self.min_time) + xmin = (datetime.datetime.combine(dt.date(), + datetime.time(int(self.xmin), 0, 0))-d1970).total_seconds() + else: + xmin = self.xmin - ''' - def deleteanotherfiles(self): - figurenames=[] - if self.figurelist != None: - for n, eachfigure in enumerate(self.figurelist): - #add specific name for each channel in channelList - ghostfigname = os.path.join(self.save, '{}_{}_{}'.format(self.titles[n].replace(' ',''),self.CODE, - datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d'))) - figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n].replace(' ',''),self.CODE, - datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S'))) - - for ghostfigure in glob.glob(ghostfigname+'*'): #ghostfigure will adopt all posible names of figures - if ghostfigure != figname: - os.remove(ghostfigure) - print 'Removing GhostFigures:' , figname - else : - '''Erasing ghost images for just on******************''' - ghostfigname = os.path.join(self.save, '{}_{}'.format(self.CODE,datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d'))) - figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S'))) - for ghostfigure in glob.glob(ghostfigname+'*'): #ghostfigure will adopt all posible names of figures - if ghostfigure != figname: - os.remove(ghostfigure) - print 'Removing GhostFigures:' , figname + if self.xmax is None: + xmax = xmin+self.xrange*60*60 + else: + if self.xaxis is 'time': + dt = datetime.datetime.fromtimestamp(self.min_time) + xmax = (datetime.datetime.combine(dt.date(), + datetime.time(int(self.xmax), 0, 0))-d1970).total_seconds() + else: + xmax = self.xmax + + ymin = self.ymin if self.ymin else numpy.nanmin(self.y) + ymax = self.ymax if self.ymax else numpy.nanmax(self.y) + + ystep = 200 if ymax>= 800 else 100 if ymax>=400 else 50 if ymax>=200 else 20 + + for n, ax in enumerate(self.axes): + if ax.firsttime: + ax.set_facecolor(self.bgcolor) + ax.yaxis.set_major_locator(MultipleLocator(ystep)) + if self.xaxis is 'time': + ax.xaxis.set_major_formatter(FuncFormatter(func)) + ax.xaxis.set_major_locator(LinearLocator(9)) + if self.xlabel is not None: + ax.set_xlabel(self.xlabel) + ax.set_ylabel(self.ylabel) + ax.firsttime = False + if self.showprofile: + self.pf_axes[n].set_ylim(ymin, ymax) + self.pf_axes[n].set_xlim(self.zmin, self.zmax) + self.pf_axes[n].set_xlabel('dB') + self.pf_axes[n].grid(b=True, axis='x') + [tick.set_visible(False) for tick in self.pf_axes[n].get_yticklabels()] + if self.colorbar: + cb = plt.colorbar(ax.plt, ax=ax, pad=0.02) + cb.ax.tick_params(labelsize=8) + if self.cb_label: + cb.set_label(self.cb_label, size=8) + elif self.cb_labels: + cb.set_label(self.cb_labels[n], size=8) + + ax.set_title('{} - {} UTC'.format( + self.titles[n], + datetime.datetime.fromtimestamp(self.max_time).strftime('%H:%M:%S')), + size=8) + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + def __plot(self): - - print 'plotting...{}'.format(self.CODE) - if self.ind_plt_ch is False : #standard + ''' + ''' + log.success('Plotting', self.name) + + self.plot() + self.format() + + for n, fig in enumerate(self.figures): + if self.nrows == 0 or self.nplots == 0: + log.warning('No data', self.name) + continue if self.show: - self.figure.show() - self.plot() - plt.tight_layout() - self.figure.canvas.manager.set_window_title('{} {} - {}'.format(self.title, self.CODE.upper(), - datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d'))) - else : - print 'len(self.figurelist): ',len(self.figurelist) - for n, eachfigure in enumerate(self.figurelist): - if self.show: - eachfigure.show() - - self.plot() - eachfigure.tight_layout() # ajuste de cada subplot - eachfigure.canvas.manager.set_window_title('{} {} - {}'.format(self.title[n], self.CODE.upper(), - datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d'))) - - # if self.save: - # if self.ind_plt_ch is False : #standard - # figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE, - # datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S'))) - # print 'Saving figure: {}'.format(figname) - # self.figure.savefig(figname) - # else : - # for n, eachfigure in enumerate(self.figurelist): - # #add specific name for each channel in channelList - # figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n],self.CODE, - # datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S'))) - # - # print 'Saving figure: {}'.format(figname) - # eachfigure.savefig(figname) - - if self.ind_plt_ch is False : - self.figure.canvas.draw() - else : - for eachfigure in self.figurelist: - eachfigure.canvas.draw() - - if self.save: - if self.ind_plt_ch is False : #standard - figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE, - datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S'))) + fig.show() + + fig.tight_layout() + fig.canvas.manager.set_window_title('{} - {}'.format(self.title, + datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d'))) + # fig.canvas.draw() + + if self.save and self.data.ended: + channels = range(self.nrows) + if self.oneFigure: + label = '' + else: + label = '_{}'.format(channels[n]) + figname = os.path.join( + self.save, + '{}{}_{}.png'.format( + self.CODE, + label, + datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S') + ) + ) print 'Saving figure: {}'.format(figname) - self.figure.savefig(figname) - else : - for n, eachfigure in enumerate(self.figurelist): - #add specific name for each channel in channelList - figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n].replace(' ',''),self.CODE, - datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S'))) - - print 'Saving figure: {}'.format(figname) - eachfigure.savefig(figname) - + fig.savefig(figname) def plot(self): - - print 'plotting...{}'.format(self.CODE.upper()) - return + ''' + ''' + raise(NotImplementedError, 'Implement this method in child class') def run(self): - print '[Starting] {}'.format(self.name) + log.success('Starting', self.name) context = zmq.Context() receiver = context.socket(zmq.SUB) @@ -212,152 +309,104 @@ class PlotData(Operation, Process): if 'server' in self.kwargs['parent']: receiver.connect('ipc:///tmp/{}.plots'.format(self.kwargs['parent']['server'])) else: - receiver.connect("ipc:///tmp/zmq.plots") - - seconds_passed = 0 + receiver.connect("ipc:///tmp/zmq.plots") while True: try: - self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)#flags=zmq.NOBLOCK - self.started = self.data['STARTED'] - self.dataOut = self.data['dataOut'] - - if (len(self.times) < len(self.data['times']) and not self.started and self.data['ENDED']): - continue - - self.times = self.data['times'] - self.times.sort() - self.throttle_value = self.data['throttle'] - self.min_time = self.times[0] - self.max_time = self.times[-1] + self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK) + + self.min_time = self.data.times[0] + self.max_time = self.data.times[-1] if self.isConfig is False: - print 'setting up' - self.setup() + self.__setup() self.isConfig = True - self.__plot() - - if self.data['ENDED'] is True: - print '********GRAPHIC ENDED********' - self.ended = True - self.isConfig = False - self.__plot() - self.deleteanotherfiles() #CLPDG - elif seconds_passed >= self.data['throttle']: - print 'passed', seconds_passed - self.__plot() - seconds_passed = 0 + + self.__plot() except zmq.Again as e: - print 'Waiting for data...' - plt.pause(2) - seconds_passed += 2 + log.log('Waiting for data...') + if self.data: + plt.pause(self.data.throttle) + else: + time.sleep(2) def close(self): - if self.dataOut: + if self.data: self.__plot() class PlotSpectraData(PlotData): + ''' + Plot for Spectra data + ''' CODE = 'spc' - colormap = 'jro' - CONFLATE = False + colormap = 'jro' def setup(self): - - ncolspan = 1 - colspan = 1 - self.ncols = int(numpy.sqrt(self.dataOut.nChannels)+0.9) - self.nrows = int(self.dataOut.nChannels*1./self.ncols + 0.9) - self.width = 3.6*self.ncols - self.height = 3.2*self.nrows - if self.showprofile: - ncolspan = 3 - colspan = 2 - self.width += 1.2*self.ncols + self.nplots = len(self.data.channels) + self.ncols = int(numpy.sqrt(self.nplots)+ 0.9) + self.nrows = int((1.0*self.nplots/self.ncols) + 0.9) + self.width = 3.4*self.ncols + self.height = 3*self.nrows + self.cb_label = 'dB' + if self.showprofile: + self.width += 0.8*self.ncols self.ylabel = 'Range [Km]' - self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList] - - if self.figure is None: - self.figure = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') - else: - self.figure.clf() - - n = 0 - for y in range(self.nrows): - for x in range(self.ncols): - if n >= self.dataOut.nChannels: - break - ax = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan), 1, colspan) - if self.showprofile: - ax.ax_profile = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan+colspan), 1, 1) - - ax.firsttime = True - self.axes.append(ax) - n += 1 def plot(self): - if self.xaxis == "frequency": - x = self.dataOut.getFreqRange(1)/1000. - xlabel = "Frequency (kHz)" + x = self.data.xrange[0] + self.xlabel = "Frequency (kHz)" elif self.xaxis == "time": - x = self.dataOut.getAcfRange(1) - xlabel = "Time (ms)" + x = self.data.xrange[1] + self.xlabel = "Time (ms)" else: - x = self.dataOut.getVelRange(1) - xlabel = "Velocity (m/s)" + x = self.data.xrange[2] + self.xlabel = "Velocity (m/s)" + + if self.CODE == 'spc_mean': + x = self.data.xrange[2] + self.xlabel = "Velocity (m/s)" - y = self.dataOut.getHeiRange() - z = self.data[self.CODE] + self.titles = [] + y = self.data.heights + self.y = y + z = self.data['spc'] + for n, ax in enumerate(self.axes): + noise = self.data['noise'][n][-1] + if self.CODE == 'spc_mean': + mean = self.data['mean'][n][-1] if ax.firsttime: - self.xmax = self.xmax if self.xmax else np.nanmax(x) + self.xmax = self.xmax if self.xmax else numpy.nanmax(x) self.xmin = self.xmin if self.xmin else -self.xmax - self.ymin = self.ymin if self.ymin else np.nanmin(y) - self.ymax = self.ymax if self.ymax else np.nanmax(y) - self.zmin = self.zmin if self.zmin else np.nanmin(z) - self.zmax = self.zmax if self.zmax else np.nanmax(z) - ax.plot = ax.pcolormesh(x, y, z[n].T, - vmin=self.zmin, - vmax=self.zmax, - cmap=plt.get_cmap(self.colormap) - ) - divider = make_axes_locatable(ax) - cax = divider.new_horizontal(size='3%', pad=0.05) - self.figure.add_axes(cax) - plt.colorbar(ax.plot, cax) - - ax.set_xlim(self.xmin, self.xmax) - ax.set_ylim(self.ymin, self.ymax) - - ax.set_ylabel(self.ylabel) - ax.set_xlabel(xlabel) - - ax.firsttime = False + self.zmin = self.zmin if self.zmin else numpy.nanmin(z) + self.zmax = self.zmax if self.zmax else numpy.nanmax(z) + ax.plt = ax.pcolormesh(x, y, z[n].T, + vmin=self.zmin, + vmax=self.zmax, + cmap=plt.get_cmap(self.colormap) + ) if self.showprofile: - ax.plot_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0] - ax.ax_profile.set_xlim(self.zmin, self.zmax) - ax.ax_profile.set_ylim(self.ymin, self.ymax) - ax.ax_profile.set_xlabel('dB') - ax.ax_profile.grid(b=True, axis='x') - ax.plot_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y, - color="k", linestyle="dashed", lw=2)[0] - [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()] + ax.plt_profile= self.pf_axes[n].plot(self.data['rti'][n][-1], y)[0] + ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y, + color="k", linestyle="dashed", lw=1)[0] + if self.CODE == 'spc_mean': + ax.plt_mean = ax.plot(mean, y, color='k')[0] else: - ax.plot.set_array(z[n].T.ravel()) + ax.plt.set_array(z[n].T.ravel()) if self.showprofile: - ax.plot_profile.set_data(self.data['rti'][self.max_time][n], y) - ax.plot_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y) + ax.plt_profile.set_data(self.data['rti'][n][-1], y) + ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y) + if self.CODE == 'spc_mean': + ax.plt_mean.set_data(mean, y) - ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]), - size=8) + self.titles.append('CH {}: {:3.2f}dB'.format(n, noise)) self.saveTime = self.max_time @@ -367,545 +416,245 @@ class PlotCrossSpectraData(PlotData): zmin_coh = None zmax_coh = None zmin_phase = None - zmax_phase = None - CONFLATE = False + zmax_phase = None def setup(self): - ncolspan = 1 - colspan = 1 - self.ncols = 2 - self.nrows = self.dataOut.nPairs - self.width = 3.6*self.ncols - self.height = 3.2*self.nrows - + self.ncols = 4 + self.nrows = len(self.data.pairs) + self.nplots = self.nrows*4 + self.width = 3.4*self.ncols + self.height = 3*self.nrows self.ylabel = 'Range [Km]' - self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList] - - if self.figure is None: - self.figure = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') - else: - self.figure.clf() - - for y in range(self.nrows): - for x in range(self.ncols): - ax = plt.subplot2grid((self.nrows, self.ncols), (y, x), 1, 1) - ax.firsttime = True - self.axes.append(ax) + self.showprofile = False def plot(self): if self.xaxis == "frequency": - x = self.dataOut.getFreqRange(1)/1000. - xlabel = "Frequency (kHz)" + x = self.data.xrange[0] + self.xlabel = "Frequency (kHz)" elif self.xaxis == "time": - x = self.dataOut.getAcfRange(1) - xlabel = "Time (ms)" + x = self.data.xrange[1] + self.xlabel = "Time (ms)" else: - x = self.dataOut.getVelRange(1) - xlabel = "Velocity (m/s)" + x = self.data.xrange[2] + self.xlabel = "Velocity (m/s)" + + self.titles = [] - y = self.dataOut.getHeiRange() - z_coh = self.data['cspc_coh'] - z_phase = self.data['cspc_phase'] + y = self.data.heights + self.y = y + spc = self.data['spc'] + cspc = self.data['cspc'] for n in range(self.nrows): - ax = self.axes[2*n] - ax1 = self.axes[2*n+1] + noise = self.data['noise'][n][-1] + pair = self.data.pairs[n] + ax = self.axes[4*n] + ax3 = self.axes[4*n+3] if ax.firsttime: - self.xmax = self.xmax if self.xmax else np.nanmax(x) + self.xmax = self.xmax if self.xmax else numpy.nanmax(x) self.xmin = self.xmin if self.xmin else -self.xmax - self.ymin = self.ymin if self.ymin else np.nanmin(y) - self.ymax = self.ymax if self.ymax else np.nanmax(y) - self.zmin_coh = self.zmin_coh if self.zmin_coh else 0.0 - self.zmax_coh = self.zmax_coh if self.zmax_coh else 1.0 - self.zmin_phase = self.zmin_phase if self.zmin_phase else -180 - self.zmax_phase = self.zmax_phase if self.zmax_phase else 180 - - ax.plot = ax.pcolormesh(x, y, z_coh[n].T, - vmin=self.zmin_coh, - vmax=self.zmax_coh, - cmap=plt.get_cmap(self.colormap_coh) - ) - divider = make_axes_locatable(ax) - cax = divider.new_horizontal(size='3%', pad=0.05) - self.figure.add_axes(cax) - plt.colorbar(ax.plot, cax) - - ax.set_xlim(self.xmin, self.xmax) - ax.set_ylim(self.ymin, self.ymax) - - ax.set_ylabel(self.ylabel) - ax.set_xlabel(xlabel) - ax.firsttime = False - - ax1.plot = ax1.pcolormesh(x, y, z_phase[n].T, - vmin=self.zmin_phase, - vmax=self.zmax_phase, - cmap=plt.get_cmap(self.colormap_phase) - ) - divider = make_axes_locatable(ax1) - cax = divider.new_horizontal(size='3%', pad=0.05) - self.figure.add_axes(cax) - plt.colorbar(ax1.plot, cax) - - ax1.set_xlim(self.xmin, self.xmax) - ax1.set_ylim(self.ymin, self.ymax) - - ax1.set_ylabel(self.ylabel) - ax1.set_xlabel(xlabel) - ax1.firsttime = False + self.zmin = self.zmin if self.zmin else numpy.nanmin(spc) + self.zmax = self.zmax if self.zmax else numpy.nanmax(spc) + ax.plt = ax.pcolormesh(x, y, spc[pair[0]].T, + vmin=self.zmin, + vmax=self.zmax, + cmap=plt.get_cmap(self.colormap) + ) else: - ax.plot.set_array(z_coh[n].T.ravel()) - ax1.plot.set_array(z_phase[n].T.ravel()) - - ax.set_title('Coherence Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8) - ax1.set_title('Phase Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8) - self.saveTime = self.max_time - + ax.plt.set_array(spc[pair[0]].T.ravel()) + self.titles.append('CH {}: {:3.2f}dB'.format(n, noise)) -class PlotSpectraMeanData(PlotSpectraData): - - CODE = 'spc_mean' - colormap = 'jet' - - def plot(self): - - if self.xaxis == "frequency": - x = self.dataOut.getFreqRange(1)/1000. - xlabel = "Frequency (kHz)" - elif self.xaxis == "time": - x = self.dataOut.getAcfRange(1) - xlabel = "Time (ms)" - else: - x = self.dataOut.getVelRange(1) - xlabel = "Velocity (m/s)" - - y = self.dataOut.getHeiRange() - z = self.data['spc'] - mean = self.data['mean'][self.max_time] - - for n, ax in enumerate(self.axes): - - if ax.firsttime: - self.xmax = self.xmax if self.xmax else np.nanmax(x) - self.xmin = self.xmin if self.xmin else -self.xmax - self.ymin = self.ymin if self.ymin else np.nanmin(y) - self.ymax = self.ymax if self.ymax else np.nanmax(y) - self.zmin = self.zmin if self.zmin else np.nanmin(z) - self.zmax = self.zmax if self.zmax else np.nanmax(z) - ax.plt = ax.pcolormesh(x, y, z[n].T, + ax = self.axes[4*n+1] + if ax.firsttime: + ax.plt = ax.pcolormesh(x, y, spc[pair[1]].T, vmin=self.zmin, vmax=self.zmax, cmap=plt.get_cmap(self.colormap) ) - ax.plt_dop = ax.plot(mean[n], y, - color='k')[0] - - divider = make_axes_locatable(ax) - cax = divider.new_horizontal(size='3%', pad=0.05) - self.figure.add_axes(cax) - plt.colorbar(ax.plt, cax) - - ax.set_xlim(self.xmin, self.xmax) - ax.set_ylim(self.ymin, self.ymax) - - ax.set_ylabel(self.ylabel) - ax.set_xlabel(xlabel) - - ax.firsttime = False - - if self.showprofile: - ax.plt_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0] - ax.ax_profile.set_xlim(self.zmin, self.zmax) - ax.ax_profile.set_ylim(self.ymin, self.ymax) - ax.ax_profile.set_xlabel('dB') - ax.ax_profile.grid(b=True, axis='x') - ax.plt_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y, - color="k", linestyle="dashed", lw=2)[0] - [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()] else: - ax.plt.set_array(z[n].T.ravel()) - ax.plt_dop.set_data(mean[n], y) - if self.showprofile: - ax.plt_profile.set_data(self.data['rti'][self.max_time][n], y) - ax.plt_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y) + ax.plt.set_array(spc[pair[1]].T.ravel()) + self.titles.append('CH {}: {:3.2f}dB'.format(n, noise)) + + out = cspc[n]/numpy.sqrt(spc[pair[0]]*spc[pair[1]]) + coh = numpy.abs(out) + phase = numpy.arctan2(out.imag, out.real)*180/numpy.pi + + ax = self.axes[4*n+2] + if ax.firsttime: + ax.plt = ax.pcolormesh(x, y, coh.T, + vmin=0, + vmax=1, + cmap=plt.get_cmap(self.colormap_coh) + ) + else: + ax.plt.set_array(coh.T.ravel()) + self.titles.append('Coherence Ch{} * Ch{}'.format(pair[0], pair[1])) - ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]), - size=8) + ax = self.axes[4*n+3] + if ax.firsttime: + ax.plt = ax.pcolormesh(x, y, phase.T, + vmin=-180, + vmax=180, + cmap=plt.get_cmap(self.colormap_phase) + ) + else: + ax.plt.set_array(phase.T.ravel()) + self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1])) + self.saveTime = self.max_time +class PlotSpectraMeanData(PlotSpectraData): + ''' + Plot for Spectra and Mean + ''' + CODE = 'spc_mean' + colormap = 'jro' + + class PlotRTIData(PlotData): + ''' + Plot for RTI data + ''' CODE = 'rti' colormap = 'jro' def setup(self): - self.ncols = 1 - self.nrows = self.dataOut.nChannels - self.width = 10 - #TODO : arreglar la altura de la figura, esta hardcodeada. - #Se arreglo, testear! - if self.ind_plt_ch: - self.height = 3.2#*self.nrows if self.nrows<6 else 12 - else: - self.height = 2.2*self.nrows if self.nrows<6 else 12 - - ''' - if self.nrows==1: - self.height += 1 - ''' + self.xaxis = 'time' + self.ncols = 1 + self.nrows = len(self.data.channels) + self.nplots = len(self.data.channels) self.ylabel = 'Range [Km]' - self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList] - - ''' - Logica: - 1) Si la variable ind_plt_ch es True, va a crear mas de 1 figura - 2) guardamos "Figures" en una lista y "axes" en otra, quizas se deberia guardar el - axis dentro de "Figures" como un diccionario. - ''' - if self.ind_plt_ch is False: #standard mode - - if self.figure is None: #solo para la priemra vez - self.figure = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') - else: - self.figure.clf() - self.axes = [] - - - for n in range(self.nrows): - ax = self.figure.add_subplot(self.nrows, self.ncols, n+1) - #ax = self.figure(n+1) - ax.firsttime = True - self.axes.append(ax) - - else : #append one figure foreach channel in channelList - if self.figurelist == None: - self.figurelist = [] - for n in range(self.nrows): - self.figure = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') - #add always one subplot - self.figurelist.append(self.figure) - - else : # cada dia nuevo limpia el axes, pero mantiene el figure - for eachfigure in self.figurelist: - eachfigure.clf() # eliminaria todas las figuras de la lista? - self.axes = [] - - for eachfigure in self.figurelist: - ax = eachfigure.add_subplot(1,1,1) #solo 1 axis por figura - #ax = self.figure(n+1) - ax.firsttime = True - #Cada figura tiene un distinto puntero - self.axes.append(ax) - #plt.close(eachfigure) - + self.cb_label = 'dB' + self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)] def plot(self): + self.x = self.data.times + self.y = self.data.heights + self.z = self.data[self.CODE] + self.z = numpy.ma.masked_invalid(self.z) - if self.ind_plt_ch is False: #standard mode - self.x = np.array(self.times) - self.y = self.dataOut.getHeiRange() - self.z = [] - - for ch in range(self.nrows): - self.z.append([self.data[self.CODE][t][ch] for t in self.times]) - - self.z = np.array(self.z) - for n, ax in enumerate(self.axes): - x, y, z = self.fill_gaps(*self.decimate()) - if self.xmin is None: - xmin = self.min_time - else: - xmin = fromtimestamp(int(self.xmin), self.min_time) - if self.xmax is None: - xmax = xmin + self.xrange*60*60 - else: - xmax = xmin + (self.xmax - self.xmin) * 60 * 60 - self.zmin = self.zmin if self.zmin else np.min(self.z) - self.zmax = self.zmax if self.zmax else np.max(self.z) - if ax.firsttime: - self.ymin = self.ymin if self.ymin else np.nanmin(self.y) - self.ymax = self.ymax if self.ymax else np.nanmax(self.y) - plot = ax.pcolormesh(x, y, z[n].T, - vmin=self.zmin, - vmax=self.zmax, - cmap=plt.get_cmap(self.colormap) - ) - divider = make_axes_locatable(ax) - cax = divider.new_horizontal(size='2%', pad=0.05) - self.figure.add_axes(cax) - plt.colorbar(plot, cax) - ax.set_ylim(self.ymin, self.ymax) - ax.xaxis.set_major_formatter(FuncFormatter(func)) - ax.xaxis.set_major_locator(LinearLocator(6)) - ax.set_ylabel(self.ylabel) - # if self.xmin is None: - # xmin = self.min_time - # else: - # xmin = (datetime.datetime.combine(self.dataOut.datatime.date(), - # datetime.time(self.xmin, 0, 0))-d1970).total_seconds() - - ax.set_xlim(xmin, xmax) - ax.firsttime = False - else: - ax.collections.remove(ax.collections[0]) - ax.set_xlim(xmin, xmax) - plot = ax.pcolormesh(x, y, z[n].T, - vmin=self.zmin, - vmax=self.zmax, - cmap=plt.get_cmap(self.colormap) - ) - ax.set_title('{} {}'.format(self.titles[n], - datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')), - size=8) - - self.saveTime = self.min_time - else : - self.x = np.array(self.times) - self.y = self.dataOut.getHeiRange() - self.z = [] - - for ch in range(self.nrows): - self.z.append([self.data[self.CODE][t][ch] for t in self.times]) - - self.z = np.array(self.z) - for n, eachfigure in enumerate(self.figurelist): #estaba ax in axes - - x, y, z = self.fill_gaps(*self.decimate()) - xmin = self.min_time - xmax = xmin+self.xrange*60*60 - self.zmin = self.zmin if self.zmin else np.min(self.z) - self.zmax = self.zmax if self.zmax else np.max(self.z) - if self.axes[n].firsttime: - self.ymin = self.ymin if self.ymin else np.nanmin(self.y) - self.ymax = self.ymax if self.ymax else np.nanmax(self.y) - plot = self.axes[n].pcolormesh(x, y, z[n].T, - vmin=self.zmin, - vmax=self.zmax, - cmap=plt.get_cmap(self.colormap) - ) - divider = make_axes_locatable(self.axes[n]) - cax = divider.new_horizontal(size='2%', pad=0.05) - eachfigure.add_axes(cax) - #self.figure2.add_axes(cax) - plt.colorbar(plot, cax) - self.axes[n].set_ylim(self.ymin, self.ymax) - - self.axes[n].xaxis.set_major_formatter(FuncFormatter(func)) - self.axes[n].xaxis.set_major_locator(LinearLocator(6)) - - self.axes[n].set_ylabel(self.ylabel) - - if self.xmin is None: - xmin = self.min_time - else: - xmin = (datetime.datetime.combine(self.dataOut.datatime.date(), - datetime.time(self.xmin, 0, 0))-d1970).total_seconds() - - self.axes[n].set_xlim(xmin, xmax) - self.axes[n].firsttime = False - else: - self.axes[n].collections.remove(self.axes[n].collections[0]) - self.axes[n].set_xlim(xmin, xmax) - plot = self.axes[n].pcolormesh(x, y, z[n].T, - vmin=self.zmin, - vmax=self.zmax, - cmap=plt.get_cmap(self.colormap) - ) - self.axes[n].set_title('{} {}'.format(self.titles[n], - datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')), - size=8) + for n, ax in enumerate(self.axes): + x, y, z = self.fill_gaps(*self.decimate()) + self.zmin = self.zmin if self.zmin else numpy.min(self.z) + self.zmax = self.zmax if self.zmax else numpy.max(self.z) + if ax.firsttime: + ax.plt = ax.pcolormesh(x, y, z[n].T, + vmin=self.zmin, + vmax=self.zmax, + cmap=plt.get_cmap(self.colormap) + ) + if self.showprofile: + ax.plot_profile= self.pf_axes[n].plot(self.data['rti'][n][-1], self.y)[0] + ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y, + color="k", linestyle="dashed", lw=1)[0] + else: + ax.collections.remove(ax.collections[0]) + ax.plt = ax.pcolormesh(x, y, z[n].T, + vmin=self.zmin, + vmax=self.zmax, + cmap=plt.get_cmap(self.colormap) + ) + if self.showprofile: + ax.plot_profile.set_data(self.data['rti'][n][-1], self.y) + ax.plot_noise.set_data(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y) - self.saveTime = self.min_time + self.saveTime = self.min_time class PlotCOHData(PlotRTIData): + ''' + Plot for Coherence data + ''' CODE = 'coh' def setup(self): - + self.xaxis = 'time' self.ncols = 1 - self.nrows = self.dataOut.nPairs - self.width = 10 - self.height = 2.2*self.nrows if self.nrows<6 else 12 - self.ind_plt_ch = False #just for coherence and phase - if self.nrows==1: - self.height += 1 - self.ylabel = 'Range [Km]' - self.titles = ['{} Ch{} * Ch{}'.format(self.CODE.upper(), x[0], x[1]) for x in self.dataOut.pairsList] - - if self.figure is None: - self.figure = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') + self.nrows = len(self.data.pairs) + self.nplots = len(self.data.pairs) + self.ylabel = 'Range [Km]' + if self.CODE == 'coh': + self.cb_label = '' + self.titles = ['Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs] else: - self.figure.clf() - self.axes = [] + self.cb_label = 'Degrees' + self.titles = ['Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs] - for n in range(self.nrows): - ax = self.figure.add_subplot(self.nrows, self.ncols, n+1) - ax.firsttime = True - self.axes.append(ax) + +class PlotPHASEData(PlotCOHData): + ''' + Plot for Phase map data + ''' + + CODE = 'phase' + colormap = 'seismic' class PlotNoiseData(PlotData): + ''' + Plot for noise + ''' + CODE = 'noise' def setup(self): - + self.xaxis = 'time' self.ncols = 1 self.nrows = 1 - self.width = 10 - self.height = 3.2 + self.nplots = 1 self.ylabel = 'Intensity [dB]' self.titles = ['Noise'] - - if self.figure is None: - self.figure = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') - else: - self.figure.clf() - self.axes = [] - - self.ax = self.figure.add_subplot(self.nrows, self.ncols, 1) - self.ax.firsttime = True + self.colorbar = False def plot(self): - x = self.times + x = self.data.times xmin = self.min_time xmax = xmin+self.xrange*60*60 - if self.ax.firsttime: - for ch in self.dataOut.channelList: - y = [self.data[self.CODE][t][ch] for t in self.times] - self.ax.plot(x, y, lw=1, label='Ch{}'.format(ch)) - self.ax.firsttime = False - self.ax.xaxis.set_major_formatter(FuncFormatter(func)) - self.ax.xaxis.set_major_locator(LinearLocator(6)) - self.ax.set_ylabel(self.ylabel) + Y = self.data[self.CODE] + + if self.axes[0].firsttime: + for ch in self.data.channels: + y = Y[ch] + self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch)) plt.legend() else: - for ch in self.dataOut.channelList: - y = [self.data[self.CODE][t][ch] for t in self.times] - self.ax.lines[ch].set_data(x, y) - - self.ax.set_xlim(xmin, xmax) - self.ax.set_ylim(min(y)-5, max(y)+5) + for ch in self.data.channels: + y = Y[ch] + self.axes[0].lines[ch].set_data(x, y) + + self.ymin = numpy.nanmin(Y) - 5 + self.ymax = numpy.nanmax(Y) + 5 self.saveTime = self.min_time -class PlotWindProfilerData(PlotRTIData): - - CODE = 'wind' - colormap = 'seismic' - - def setup(self): - self.ncols = 1 - self.nrows = self.dataOut.data_output.shape[0] - self.width = 10 - self.height = 2.2*self.nrows - self.ylabel = 'Height [Km]' - self.titles = ['Zonal Wind' ,'Meridional Wind', 'Vertical Wind'] - self.clabels = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)'] - self.windFactor = [1, 1, 100] - - if self.figure is None: - self.figure = plt.figure(figsize=(self.width, self.height), - edgecolor='k', - facecolor='w') - else: - self.figure.clf() - self.axes = [] - - for n in range(self.nrows): - ax = self.figure.add_subplot(self.nrows, self.ncols, n+1) - ax.firsttime = True - self.axes.append(ax) - - def plot(self): - - self.x = np.array(self.times) - self.y = self.dataOut.heightList - self.z = [] - - for ch in range(self.nrows): - self.z.append([self.data['output'][t][ch] for t in self.times]) - - self.z = np.array(self.z) - self.z = numpy.ma.masked_invalid(self.z) - - cmap=plt.get_cmap(self.colormap) - cmap.set_bad('black', 1.) - - for n, ax in enumerate(self.axes): - x, y, z = self.fill_gaps(*self.decimate()) - xmin = self.min_time - xmax = xmin+self.xrange*60*60 - if ax.firsttime: - self.ymin = self.ymin if self.ymin else np.nanmin(self.y) - self.ymax = self.ymax if self.ymax else np.nanmax(self.y) - self.zmax = self.zmax if self.zmax else numpy.nanmax(abs(self.z[:-1, :])) - self.zmin = self.zmin if self.zmin else -self.zmax - - plot = ax.pcolormesh(x, y, z[n].T*self.windFactor[n], - vmin=self.zmin, - vmax=self.zmax, - cmap=cmap - ) - divider = make_axes_locatable(ax) - cax = divider.new_horizontal(size='2%', pad=0.05) - self.figure.add_axes(cax) - cb = plt.colorbar(plot, cax) - cb.set_label(self.clabels[n]) - ax.set_ylim(self.ymin, self.ymax) - - ax.xaxis.set_major_formatter(FuncFormatter(func)) - ax.xaxis.set_major_locator(LinearLocator(6)) - - ax.set_ylabel(self.ylabel) - - ax.set_xlim(xmin, xmax) - ax.firsttime = False - else: - ax.collections.remove(ax.collections[0]) - ax.set_xlim(xmin, xmax) - plot = ax.pcolormesh(x, y, z[n].T*self.windFactor[n], - vmin=self.zmin, - vmax=self.zmax, - cmap=plt.get_cmap(self.colormap) - ) - ax.set_title('{} {}'.format(self.titles[n], - datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')), - size=8) - - self.saveTime = self.min_time - - class PlotSNRData(PlotRTIData): + ''' + Plot for SNR Data + ''' + CODE = 'snr' colormap = 'jet' + class PlotDOPData(PlotRTIData): + ''' + Plot for DOPPLER Data + ''' + CODE = 'dop' colormap = 'jet' -class PlotPHASEData(PlotCOHData): - CODE = 'phase' - colormap = 'seismic' - - class PlotSkyMapData(PlotData): + ''' + Plot for meteors detection data + ''' CODE = 'met' @@ -932,7 +681,7 @@ class PlotSkyMapData(PlotData): def plot(self): - arrayParameters = np.concatenate([self.data['param'][t] for t in self.times]) + arrayParameters = numpy.concatenate([self.data['param'][t] for t in self.data.times]) error = arrayParameters[:,-1] indValid = numpy.where(error == 0)[0] finalMeteor = arrayParameters[indValid,:] @@ -962,3 +711,72 @@ class PlotSkyMapData(PlotData): self.ax.set_title(title, size=8) self.saveTime = self.max_time + +class PlotParamData(PlotRTIData): + ''' + Plot for data_param object + ''' + + CODE = 'param' + colormap = 'seismic' + + def setup(self): + self.xaxis = 'time' + self.ncols = 1 + self.nrows = self.data.shape(self.CODE)[0] + self.nplots = self.nrows + if self.showSNR: + self.nrows += 1 + + self.ylabel = 'Height [Km]' + self.titles = self.data.parameters \ + if self.data.parameters else ['Param {}'.format(x) for x in xrange(self.nrows)] + if self.showSNR: + self.titles.append('SNR') + + def plot(self): + self.data.normalize_heights() + self.x = self.data.times + self.y = self.data.heights + if self.showSNR: + self.z = numpy.concatenate( + (self.data[self.CODE], self.data['snr']) + ) + else: + self.z = self.data[self.CODE] + + self.z = numpy.ma.masked_invalid(self.z) + + for n, ax in enumerate(self.axes): + + x, y, z = self.fill_gaps(*self.decimate()) + + if ax.firsttime: + if self.zlimits is not None: + self.zmin, self.zmax = self.zlimits[n] + self.zmax = self.zmax if self.zmax is not None else numpy.nanmax(abs(self.z[:-1, :])) + self.zmin = self.zmin if self.zmin is not None else -self.zmax + ax.plt = ax.pcolormesh(x, y, z[n, :, :].T*self.factors[n], + vmin=self.zmin, + vmax=self.zmax, + cmap=self.cmaps[n] + ) + else: + if self.zlimits is not None: + self.zmin, self.zmax = self.zlimits[n] + ax.collections.remove(ax.collections[0]) + ax.plt = ax.pcolormesh(x, y, z[n, :, :].T*self.factors[n], + vmin=self.zmin, + vmax=self.zmax, + cmap=self.cmaps[n] + ) + + self.saveTime = self.min_time + +class PlotOuputData(PlotParamData): + ''' + Plot data_output object + ''' + + CODE = 'output' + colormap = 'seismic' \ No newline at end of file diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 18f2ced..8fe1cf2 100644 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -1,39 +1,79 @@ import numpy import math from scipy import optimize, interpolate, signal, stats, ndimage +import scipy import re import datetime import copy import sys import importlib import itertools - +from multiprocessing import Pool, TimeoutError +from multiprocessing.pool import ThreadPool +import copy_reg +import cPickle +import types +from functools import partial +import time +#from sklearn.cluster import KMeans + +import matplotlib.pyplot as plt + +from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters from jroproc_base import ProcessingUnit, Operation from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon +from scipy import asarray as ar,exp +from scipy.optimize import curve_fit +import warnings +from numpy import NaN +from scipy.optimize.optimize import OptimizeWarning +warnings.filterwarnings('ignore') -class ParametersProc(ProcessingUnit): +SPEED_OF_LIGHT = 299792458 + + +'''solving pickling issue''' + +def _pickle_method(method): + func_name = method.im_func.__name__ + obj = method.im_self + cls = method.im_class + return _unpickle_method, (func_name, obj, cls) + +def _unpickle_method(func_name, obj, cls): + for cls in cls.mro(): + try: + func = cls.__dict__[func_name] + except KeyError: + pass + else: + break + return func.__get__(obj, cls) + +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 @@ -55,25 +95,25 @@ class ParametersProc(ProcessingUnit): self.dataOut.ippSeconds = self.dataIn.ippSeconds # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter self.dataOut.timeInterval1 = self.dataIn.timeInterval - self.dataOut.heightList = self.dataIn.getHeiRange() + self.dataOut.heightList = self.dataIn.getHeiRange() self.dataOut.frequency = self.dataIn.frequency - #self.dataOut.noise = self.dataIn.noise - + # self.dataOut.noise = self.dataIn.noise + def run(self): - + #---------------------- Voltage Data --------------------------- - + if self.dataIn.type == "Voltage": self.__updateObjFromInput() self.dataOut.data_pre = self.dataIn.data.copy() self.dataOut.flagNoData = False self.dataOut.utctimeInit = self.dataIn.utctime - self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds + self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds return - + #---------------------- Spectra Data --------------------------- - + if self.dataIn.type == "Spectra": self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc) @@ -83,107 +123,1307 @@ class ParametersProc(ProcessingUnit): self.dataOut.nIncohInt = self.dataIn.nIncohInt self.dataOut.nFFTPoints = self.dataIn.nFFTPoints self.dataOut.ippFactor = self.dataIn.ippFactor - #self.dataOut.normFactor = self.dataIn.getNormFactor() - self.dataOut.pairsList = self.dataIn.pairsList - self.dataOut.groupList = self.dataIn.pairsList self.dataOut.abscissaList = self.dataIn.getVelRange(1) + self.dataOut.spc_noise = self.dataIn.getNoise() + self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1)) + self.dataOut.pairsList = self.dataIn.pairsList + self.dataOut.groupList = self.dataIn.pairsList self.dataOut.flagNoData = False - + + if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels + self.dataOut.ChanDist = self.dataIn.ChanDist + else: self.dataOut.ChanDist = None + + if hasattr(self.dataIn, 'VelRange'): #Velocities range + self.dataOut.VelRange = self.dataIn.VelRange + else: self.dataOut.VelRange = None + + if hasattr(self.dataIn, 'RadarConst'): #Radar Constant + self.dataOut.RadarConst = self.dataIn.RadarConst + + if hasattr(self.dataIn, 'NPW'): #NPW + self.dataOut.NPW = self.dataIn.NPW + + if hasattr(self.dataIn, 'COFA'): #COFA + self.dataOut.COFA = self.dataIn.COFA + + + #---------------------- Correlation Data --------------------------- - + if self.dataIn.type == "Correlation": acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions() - + self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:]) self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:]) self.dataOut.groupList = (acf_pairs, ccf_pairs) - + self.dataOut.abscissaList = self.dataIn.lagRange self.dataOut.noise = self.dataIn.noise self.dataOut.data_SNR = self.dataIn.SNR self.dataOut.flagNoData = False self.dataOut.nAvg = self.dataIn.nAvg - + #---------------------- Parameters Data --------------------------- - + if self.dataIn.type == "Parameters": self.dataOut.copy(self.dataIn) - self.dataOut.utctimeInit = self.dataIn.utctime self.dataOut.flagNoData = False - + return True - + self.__updateObjFromInput() self.dataOut.utctimeInit = self.dataIn.utctime self.dataOut.paramInterval = self.dataIn.timeInterval - + return -class SpectralMoments(Operation): +def target(tups): + + obj, args = tups + #print 'TARGETTT', obj, args + return obj.FitGau(args) + +class GaussianFit(Operation): + ''' - Function SpectralMoments() + Function that fit of one and two generalized gaussians (gg) based + on the PSD shape across an "power band" identified from a cumsum of + the measured spectrum - noise. + + Input: + self.dataOut.data_pre : SelfSpectra + + Output: + self.dataOut.GauSPC : SPC_ch1, SPC_ch2 + + ''' + def __init__(self, **kwargs): + Operation.__init__(self, **kwargs) + self.i=0 + + + def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points + """This routine will find a couple of generalized Gaussians to a power spectrum + input: spc + output: + Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise + """ + + self.spc = dataOut.data_pre[0].copy() + + + print 'SelfSpectra Shape', numpy.asarray(self.spc).shape + + + #plt.figure(50) + #plt.subplot(121) + #plt.plot(self.spc,'k',label='spc(66)') + #plt.plot(xFrec,ySamples[1],'g',label='Ch1') + #plt.plot(xFrec,ySamples[2],'r',label='Ch2') + #plt.plot(xFrec,FitGauss,'yo:',label='fit') + #plt.legend() + #plt.title('DATOS A ALTURA DE 7500 METROS') + #plt.show() + + self.Num_Hei = self.spc.shape[2] + #self.Num_Bin = len(self.spc) + self.Num_Bin = self.spc.shape[1] + self.Num_Chn = self.spc.shape[0] + + Vrange = dataOut.abscissaList + + #print 'self.spc2', numpy.asarray(self.spc).shape + + GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei]) + SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) + SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) + SPC_ch1[:] = numpy.NaN + SPC_ch2[:] = numpy.NaN - Calculates moments (power, mean, standard deviation) and SNR of the signal + + start_time = time.time() + + noise_ = dataOut.spc_noise[0].copy() + + + + pool = Pool(processes=self.Num_Chn) + args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)] + objs = [self for __ in range(self.Num_Chn)] + attrs = zip(objs, args) + gauSPC = pool.map(target, attrs) + dataOut.GauSPC = numpy.asarray(gauSPC) +# ret = [] +# for n in range(self.Num_Chn): +# self.FitGau(args[n]) +# dataOut.GauSPC = ret + + + +# for ch in range(self.Num_Chn): +# +# for ht in range(self.Num_Hei): +# #print (numpy.asarray(self.spc).shape) +# spc = numpy.asarray(self.spc)[ch,:,ht] +# +# ############################################# +# # normalizing spc and noise +# # This part differs from gg1 +# spc_norm_max = max(spc) +# spc = spc / spc_norm_max +# pnoise = pnoise / spc_norm_max +# ############################################# +# +# if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's +# fatspectra=1.0 +# else: +# fatspectra=0.5 +# +# wnoise = noise_ / spc_norm_max +# #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise +# #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used +# #if wnoise>1.1*pnoise: # to be tested later +# # wnoise=pnoise +# noisebl=wnoise*0.9; noisebh=wnoise*1.1 +# spc=spc-wnoise +# +# minx=numpy.argmin(spc) +# spcs=numpy.roll(spc,-minx) +# cum=numpy.cumsum(spcs) +# tot_noise=wnoise * self.Num_Bin #64; +# #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? ''' +# #snr=tot_signal/tot_noise +# #snr=cum[-1]/tot_noise +# +# #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise +# +# snr = sum(spcs)/tot_noise +# snrdB=10.*numpy.log10(snr) +# +# #if snrdB < -9 : +# # snrdB = numpy.NaN +# # continue +# +# #print 'snr',snrdB # , sum(spcs) , tot_noise +# +# +# #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: +# # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None +# +# cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region +# cumlo=cummax*epsi; +# cumhi=cummax*(1-epsi) +# powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum-9: # when SNR is strong pick the peak with least shift (LOS velocity) error +# if oneG: +# choice=0 +# else: +# w1=lsq2[0][1]; w2=lsq2[0][5] +# a1=lsq2[0][2]; a2=lsq2[0][6] +# p1=lsq2[0][3]; p2=lsq2[0][7] +# s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2; +# gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling +# +# if gp1>gp2: +# if a1>0.7*a2: +# choice=1 +# else: +# choice=2 +# elif gp2>gp1: +# if a2>0.7*a1: +# choice=2 +# else: +# choice=1 +# else: +# choice=numpy.argmax([a1,a2])+1 +# #else: +# #choice=argmin([std2a,std2b])+1 +# +# else: # with low SNR go to the most energetic peak +# choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]]) +# +# #print 'choice',choice +# +# if choice==0: # pick the single gaussian fit +# Amplitude0=lsq1[0][2] +# shift0=lsq1[0][0] +# width0=lsq1[0][1] +# p0=lsq1[0][3] +# Amplitude1=0. +# shift1=0. +# width1=0. +# p1=0. +# noise=lsq1[0][4] +# elif choice==1: # take the first one of the 2 gaussians fitted +# Amplitude0 = lsq2[0][2] +# shift0 = lsq2[0][0] +# width0 = lsq2[0][1] +# p0 = lsq2[0][3] +# Amplitude1 = lsq2[0][6] # This is 0 in gg1 +# shift1 = lsq2[0][4] # This is 0 in gg1 +# width1 = lsq2[0][5] # This is 0 in gg1 +# p1 = lsq2[0][7] # This is 0 in gg1 +# noise = lsq2[0][8] +# else: # the second one +# Amplitude0 = lsq2[0][6] +# shift0 = lsq2[0][4] +# width0 = lsq2[0][5] +# p0 = lsq2[0][7] +# Amplitude1 = lsq2[0][2] # This is 0 in gg1 +# shift1 = lsq2[0][0] # This is 0 in gg1 +# width1 = lsq2[0][1] # This is 0 in gg1 +# p1 = lsq2[0][3] # This is 0 in gg1 +# noise = lsq2[0][8] +# +# #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0) +# SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0 +# SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1 +# #print 'SPC_ch1.shape',SPC_ch1.shape +# #print 'SPC_ch2.shape',SPC_ch2.shape +# #dataOut.data_param = SPC_ch1 +# GauSPC[0] = SPC_ch1 +# GauSPC[1] = SPC_ch2 + +# #plt.gcf().clear() +# plt.figure(50+self.i) +# self.i=self.i+1 +# #plt.subplot(121) +# plt.plot(self.spc,'k')#,label='spc(66)') +# plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1') +# #plt.plot(SPC_ch2,'r')#,label='gg2') +# #plt.plot(xFrec,ySamples[1],'g',label='Ch1') +# #plt.plot(xFrec,ySamples[2],'r',label='Ch2') +# #plt.plot(xFrec,FitGauss,'yo:',label='fit') +# plt.legend() +# plt.title('DATOS A ALTURA DE 7500 METROS') +# plt.show() +# print 'shift0', shift0 +# print 'Amplitude0', Amplitude0 +# print 'width0', width0 +# print 'p0', p0 +# print '========================' +# print 'shift1', shift1 +# print 'Amplitude1', Amplitude1 +# print 'width1', width1 +# print 'p1', p1 +# print 'noise', noise +# print 's_noise', wnoise + + print '========================================================' + print 'total_time: ', time.time()-start_time + + # re-normalizing spc and noise + # This part differs from gg1 + + + + ''' Parameters: + 1. Amplitude + 2. Shift + 3. Width + 4. Power + ''' + + + ############################################################################### + def FitGau(self, X): + + Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X + #print 'VARSSSS', ch, pnoise, noise, num_intg + + #print 'HEIGHTS', self.Num_Hei + + GauSPC = [] + SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) + SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) + SPC_ch1[:] = 0#numpy.NaN + SPC_ch2[:] = 0#numpy.NaN + + + + for ht in range(self.Num_Hei): + #print (numpy.asarray(self.spc).shape) + + #print 'TTTTT', ch , ht + #print self.spc.shape + + + spc = numpy.asarray(self.spc)[ch,:,ht] + + ############################################# + # normalizing spc and noise + # This part differs from gg1 + spc_norm_max = max(spc) + spc = spc / spc_norm_max + pnoise = pnoise / spc_norm_max + ############################################# + + fatspectra=1.0 + + wnoise = noise_ / spc_norm_max + #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used + #if wnoise>1.1*pnoise: # to be tested later + # wnoise=pnoise + noisebl=wnoise*0.9; noisebh=wnoise*1.1 + spc=spc-wnoise + # print 'wnoise', noise_[0], spc_norm_max, wnoise + minx=numpy.argmin(spc) + spcs=numpy.roll(spc,-minx) + cum=numpy.cumsum(spcs) + tot_noise=wnoise * self.Num_Bin #64; + #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise + #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? ''' + #snr=tot_signal/tot_noise + #snr=cum[-1]/tot_noise + snr = sum(spcs)/tot_noise + snrdB=10.*numpy.log10(snr) + + if snrdB < SNRlimit : + snr = numpy.NaN + SPC_ch1[:,ht] = 0#numpy.NaN + SPC_ch1[:,ht] = 0#numpy.NaN + GauSPC = (SPC_ch1,SPC_ch2) + continue + #print 'snr',snrdB #, sum(spcs) , tot_noise + + + + #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: + # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None + + cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region + cumlo=cummax*epsi; + cumhi=cummax*(1-epsi) + powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum-6: # when SNR is strong pick the peak with least shift (LOS velocity) error + if oneG: + choice=0 + else: + w1=lsq2[0][1]; w2=lsq2[0][5] + a1=lsq2[0][2]; a2=lsq2[0][6] + p1=lsq2[0][3]; p2=lsq2[0][7] + s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; + s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2; + gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling + + if gp1>gp2: + if a1>0.7*a2: + choice=1 + else: + choice=2 + elif gp2>gp1: + if a2>0.7*a1: + choice=2 + else: + choice=1 + else: + choice=numpy.argmax([a1,a2])+1 + #else: + #choice=argmin([std2a,std2b])+1 + + else: # with low SNR go to the most energetic peak + choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]]) + + + shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0]) + shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0]) + + max_vel = 20 + + #first peak will be 0, second peak will be 1 + if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range + shift0=lsq2[0][0] + width0=lsq2[0][1] + Amplitude0=lsq2[0][2] + p0=lsq2[0][3] + + shift1=lsq2[0][4] + width1=lsq2[0][5] + Amplitude1=lsq2[0][6] + p1=lsq2[0][7] + noise=lsq2[0][8] + else: + shift1=lsq2[0][0] + width1=lsq2[0][1] + Amplitude1=lsq2[0][2] + p1=lsq2[0][3] + + shift0=lsq2[0][4] + width0=lsq2[0][5] + Amplitude0=lsq2[0][6] + p0=lsq2[0][7] + noise=lsq2[0][8] + + if Amplitude0<0.1: # in case the peak is noise + shift0,width0,Amplitude0,p0 = 4*[numpy.NaN] + if Amplitude1<0.1: + shift1,width1,Amplitude1,p1 = 4*[numpy.NaN] + + +# if choice==0: # pick the single gaussian fit +# Amplitude0=lsq1[0][2] +# shift0=lsq1[0][0] +# width0=lsq1[0][1] +# p0=lsq1[0][3] +# Amplitude1=0. +# shift1=0. +# width1=0. +# p1=0. +# noise=lsq1[0][4] +# elif choice==1: # take the first one of the 2 gaussians fitted +# Amplitude0 = lsq2[0][2] +# shift0 = lsq2[0][0] +# width0 = lsq2[0][1] +# p0 = lsq2[0][3] +# Amplitude1 = lsq2[0][6] # This is 0 in gg1 +# shift1 = lsq2[0][4] # This is 0 in gg1 +# width1 = lsq2[0][5] # This is 0 in gg1 +# p1 = lsq2[0][7] # This is 0 in gg1 +# noise = lsq2[0][8] +# else: # the second one +# Amplitude0 = lsq2[0][6] +# shift0 = lsq2[0][4] +# width0 = lsq2[0][5] +# p0 = lsq2[0][7] +# Amplitude1 = lsq2[0][2] # This is 0 in gg1 +# shift1 = lsq2[0][0] # This is 0 in gg1 +# width1 = lsq2[0][1] # This is 0 in gg1 +# p1 = lsq2[0][3] # This is 0 in gg1 +# noise = lsq2[0][8] + + #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0) + SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0 + SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1 + #print 'SPC_ch1.shape',SPC_ch1.shape + #print 'SPC_ch2.shape',SPC_ch2.shape + #dataOut.data_param = SPC_ch1 + GauSPC = (SPC_ch1,SPC_ch2) + #GauSPC[1] = SPC_ch2 + +# print 'shift0', shift0 +# print 'Amplitude0', Amplitude0 +# print 'width0', width0 +# print 'p0', p0 +# print '========================' +# print 'shift1', shift1 +# print 'Amplitude1', Amplitude1 +# print 'width1', width1 +# print 'p1', p1 +# print 'noise', noise +# print 's_noise', wnoise + + return GauSPC + + + def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination. + y_model=self.y_model1(x,state) + s0,w0,a0,p0,n=state + e0=((x-s0)/w0)**2; + + e0u=((x-s0-self.Num_Bin)/w0)**2; + + e0d=((x-s0+self.Num_Bin)/w0)**2 + m0=numpy.exp(-0.5*e0**(p0/2.)); + m0u=numpy.exp(-0.5*e0u**(p0/2.)); + m0d=numpy.exp(-0.5*e0d**(p0/2.)) + JA=m0+m0u+m0d + JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d) + + JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0) + + JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2 + jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model]) + return jack1.T + + def y_jacobian2(self,x,state): + y_model=self.y_model2(x,state) + s0,w0,a0,p0,s1,w1,a1,p1,n=state + e0=((x-s0)/w0)**2; + + e0u=((x-s0- self.Num_Bin )/w0)**2; + + e0d=((x-s0+ self.Num_Bin )/w0)**2 + e1=((x-s1)/w1)**2; + + e1u=((x-s1- self.Num_Bin )/w1)**2; + + e1d=((x-s1+ self.Num_Bin )/w1)**2 + m0=numpy.exp(-0.5*e0**(p0/2.)); + m0u=numpy.exp(-0.5*e0u**(p0/2.)); + m0d=numpy.exp(-0.5*e0d**(p0/2.)) + m1=numpy.exp(-0.5*e1**(p1/2.)); + m1u=numpy.exp(-0.5*e1u**(p1/2.)); + m1d=numpy.exp(-0.5*e1d**(p1/2.)) + JA=m0+m0u+m0d + JA1=m1+m1u+m1d + JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d) + JP1=(-1/4.)*a1*m1*e1**(p1/2.)*numpy.log(e1)+(-1/4.)*a1*m1u*e1u**(p1/2.)*numpy.log(e1u)+(-1/4.)*a1*m1d*e1d**(p1/2.)*numpy.log(e1d) + + JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0) + + JS1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1) + + JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2 + + JW1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)**2+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)**2+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)**2 + jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model]) + return jack2.T + + def y_model1(self,x,state): + shift0,width0,amplitude0,power0,noise=state + model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0) + + model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0) + + model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0) + return model0+model0u+model0d+noise + + def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist + shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state + model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0) + + model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0) + + model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0) + model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1) + + model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1) + + model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1) + return model0+model0u+model0d+model1+model1u+model1d+noise + + def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is. - Type of dataIn: Spectra + return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))**2)#/(64-5.) # /(64-5.) can be commented + + def misfit2(self,state,y_data,x,num_intg): + return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.) + - Configuration Parameters: +class PrecipitationProc(Operation): + + ''' + Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R) + + Input: + self.dataOut.data_pre : SelfSpectra + + Output: + + self.dataOut.data_output : Reflectivity factor, rainfall Rate + + + Parameters affected: + ''' + + + def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None, + tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None): + + self.spc = dataOut.data_pre[0].copy() + self.Num_Hei = self.spc.shape[2] + self.Num_Bin = self.spc.shape[1] + self.Num_Chn = self.spc.shape[0] + + Velrange = dataOut.abscissaList + + if radar == "MIRA35C" : + + Ze = self.dBZeMODE2(dataOut) + + else: + + self.Pt = Pt + self.Gt = Gt + self.Gr = Gr + self.Lambda = Lambda + self.aL = aL + self.tauW = tauW + self.ThetaT = ThetaT + self.ThetaR = ThetaR + + RadarConstant = GetRadarConstant() + SPCmean = numpy.mean(self.spc,0) + ETA = numpy.zeros(self.Num_Hei) + Pr = numpy.sum(SPCmean,0) + + #for R in range(self.Num_Hei): + # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA) + + D_range = numpy.zeros(self.Num_Hei) + EqSec = numpy.zeros(self.Num_Hei) + del_V = numpy.zeros(self.Num_Hei) + + for R in range(self.Num_Hei): + ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA) + + h = R + Altitude #Range from ground to radar pulse altitude + del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity + + D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity + SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma) + + N_dist[R] = ETA[R] / SIGMA[R] + + Ze = (ETA * Lambda**4) / (numpy.pi * Km) + Z = numpy.sum( N_dist * D_range**6 ) + RR = 6*10**-4*numpy.pi * numpy.sum( D_range**3 * N_dist * Velrange ) #Rainfall rate + + + RR = (Ze/200)**(1/1.6) + dBRR = 10*numpy.log10(RR) + + dBZe = 10*numpy.log10(Ze) + dataOut.data_output = Ze + dataOut.data_param = numpy.ones([2,self.Num_Hei]) + dataOut.channelList = [0,1] + print 'channelList', dataOut.channelList + dataOut.data_param[0]=dBZe + dataOut.data_param[1]=dBRR + print 'RR SHAPE', dBRR.shape + print 'Ze SHAPE', dBZe.shape + print 'dataOut.data_param SHAPE', dataOut.data_param.shape + + + def dBZeMODE2(self, dataOut): # Processing for MIRA35C + + NPW = dataOut.NPW + COFA = dataOut.COFA + + SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]]) + RadarConst = dataOut.RadarConst + #frequency = 34.85*10**9 + + ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei])) + data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN + + ETA = numpy.sum(SNR,1) + print 'ETA' , ETA + ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN) + + Ze = numpy.ones([self.Num_Chn, self.Num_Hei] ) + + for r in range(self.Num_Hei): + + Ze[0,r] = ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2) + #Ze[1,r] = ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2) + + return Ze + + def GetRadarConstant(self): + + """ + Constants: + + Pt: Transmission Power dB + Gt: Transmission Gain dB + Gr: Reception Gain dB + Lambda: Wavelenght m + aL: Attenuation loses dB + tauW: Width of transmission pulse s + ThetaT: Transmission antenna bean angle rad + ThetaR: Reception antenna beam angle rad + + """ + Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) + Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR) + RadarConstant = Numerator / Denominator + + return RadarConstant + + + +class FullSpectralAnalysis(Operation): + + """ + Function that implements Full Spectral Analisys technique. + + Input: + self.dataOut.data_pre : SelfSpectra and CrossSPectra data + self.dataOut.groupList : Pairlist of channels + self.dataOut.ChanDist : Physical distance between receivers + + + Output: + + self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind + + + Parameters affected: Winds, height range, SNR + + """ + def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7): + + spc = dataOut.data_pre[0].copy() + cspc = dataOut.data_pre[1].copy() + + nChannel = spc.shape[0] + nProfiles = spc.shape[1] + nHeights = spc.shape[2] + + pairsList = dataOut.groupList + if dataOut.ChanDist is not None : + ChanDist = dataOut.ChanDist + else: + ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]]) + + #print 'ChanDist', ChanDist + + if dataOut.VelRange is not None: + VelRange= dataOut.VelRange + else: + VelRange= dataOut.abscissaList + + ySamples=numpy.ones([nChannel,nProfiles]) + phase=numpy.ones([nChannel,nProfiles]) + CSPCSamples=numpy.ones([nChannel,nProfiles],dtype=numpy.complex_) + coherence=numpy.ones([nChannel,nProfiles]) + PhaseSlope=numpy.ones(nChannel) + PhaseInter=numpy.ones(nChannel) + dataSNR = dataOut.data_SNR + + + + data = dataOut.data_pre + noise = dataOut.noise + print 'noise',noise + #SNRdB = 10*numpy.log10(dataOut.data_SNR) + + FirstMoment = numpy.average(dataOut.data_param[:,1,:],0) + #SNRdBMean = [] + + #for j in range(nHeights): + # FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]])) + # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]])) + + data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN + + velocityX=[] + velocityY=[] + velocityV=[] + + dbSNR = 10*numpy.log10(dataSNR) + dbSNR = numpy.average(dbSNR,0) + for Height in range(nHeights): + + [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit) + + if abs(Vzon)<100. and abs(Vzon)> 0.: + velocityX=numpy.append(velocityX, Vzon)#Vmag + + else: + print 'Vzon',Vzon + velocityX=numpy.append(velocityX, numpy.NaN) + + if abs(Vmer)<100. and abs(Vmer) > 0.: + velocityY=numpy.append(velocityY, Vmer)#Vang + + else: + print 'Vmer',Vmer + velocityY=numpy.append(velocityY, numpy.NaN) + + if dbSNR[Height] > SNRlimit: + velocityV=numpy.append(velocityV, FirstMoment[Height]) + else: + velocityV=numpy.append(velocityV, numpy.NaN) + #FirstMoment[Height]= numpy.NaN +# if SNRdBMean[Height] <12: +# FirstMoment[Height] = numpy.NaN +# velocityX[Height] = numpy.NaN +# velocityY[Height] = numpy.NaN + + + data_output[0]=numpy.array(velocityX) + data_output[1]=numpy.array(velocityY) + data_output[2]=-velocityV#FirstMoment + + print ' ' + #print 'FirstMoment' + #print FirstMoment + print 'velocityX',data_output[0] + print ' ' + print 'velocityY',data_output[1] + #print numpy.array(velocityY) + print ' ' + #print 'SNR' + #print 10*numpy.log10(dataOut.data_SNR) + #print numpy.shape(10*numpy.log10(dataOut.data_SNR)) + print ' ' + + + dataOut.data_output=data_output + return + + + def moving_average(self,x, N=2): + return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):] + + def gaus(self,xSamples,a,x0,sigma): + return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2)) + + def Find(self,x,value): + for index in range(len(x)): + if x[index]==value: + return index + + def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit): + + ySamples=numpy.ones([spc.shape[0],spc.shape[1]]) + phase=numpy.ones([spc.shape[0],spc.shape[1]]) + CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) + coherence=numpy.ones([spc.shape[0],spc.shape[1]]) + PhaseSlope=numpy.ones(spc.shape[0]) + PhaseInter=numpy.ones(spc.shape[0]) + xFrec=VelRange + + '''Getting Eij and Nij''' + + E01=ChanDist[0][0] + N01=ChanDist[0][1] + + E02=ChanDist[1][0] + N02=ChanDist[1][1] + + E12=ChanDist[2][0] + N12=ChanDist[2][1] + + z = spc.copy() + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + + for i in range(spc.shape[0]): + + '''****** Line of Data SPC ******''' + zline=z[i,:,Height] + + '''****** SPC is normalized ******''' + FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy()) + FactNorm= FactNorm/numpy.sum(FactNorm) + + SmoothSPC=self.moving_average(FactNorm,N=3) + + xSamples = ar(range(len(SmoothSPC))) + ySamples[i] = SmoothSPC + + #dbSNR=10*numpy.log10(dataSNR) + print ' ' + print ' ' + print ' ' + + #print 'dataSNR', dbSNR.shape, dbSNR[0,40:120] + print 'SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20] + print 'noise',noise + print 'zline',zline.shape, zline[0:20] + print 'FactNorm',FactNorm.shape, FactNorm[0:20] + print 'FactNorm suma', numpy.sum(FactNorm) + + for i in range(spc.shape[0]): + + '''****** Line of Data CSPC ******''' + cspcLine=cspc[i,:,Height].copy() + + '''****** CSPC is normalized ******''' + chan_index0 = pairsList[i][0] + chan_index1 = pairsList[i][1] + CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) # + + CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor) + + CSPCSamples[i] = CSPCNorm + coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor) + + coherence[i]= self.moving_average(coherence[i],N=2) + + phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi + + print 'cspcLine', cspcLine.shape, cspcLine[0:20] + print 'CSPCFactor', CSPCFactor#, CSPCFactor[0:20] + print numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i] + print 'CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20] + print 'CSPCNorm suma', numpy.sum(CSPCNorm) + print 'CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20] + + '''****** Getting fij width ******''' + + yMean=[] + yMean2=[] + + for j in range(len(ySamples[1])): + yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]])) + + '''******* Getting fitting Gaussian ******''' + meanGauss=sum(xSamples*yMean) / len(xSamples) + sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) + + print '****************************' + print 'len(xSamples): ',len(xSamples) + print 'yMean: ', yMean.shape, yMean[0:20] + print 'ySamples', ySamples.shape, ySamples[0,0:20] + print 'xSamples: ',xSamples.shape, xSamples[0:20] + + print 'meanGauss',meanGauss + print 'sigma',sigma + + #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001): + if dbSNR > SNRlimit : + try: + popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) + + if numpy.amax(popt)>numpy.amax(yMean)*0.3: + FitGauss=self.gaus(xSamples,*popt) + + else: + FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) + print 'Verificador: Dentro', Height + except :#RuntimeError: + FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) + + + else: + FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) + + Maximun=numpy.amax(yMean) + eMinus1=Maximun*numpy.exp(-1)#*0.8 + + HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1))) + HalfWidth= xFrec[HWpos] + GCpos=self.Find(FitGauss, numpy.amax(FitGauss)) + Vpos=self.Find(FactNorm, numpy.amax(FactNorm)) + + #Vpos=FirstMoment[] + + '''****** Getting Fij ******''' + + GaussCenter=xFrec[GCpos] + if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): + Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 + else: + Fij=abs(GaussCenter-HalfWidth)+0.0000001 + + '''****** Getting Frecuency range of significant data ******''' + + Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10))) + + if Rangpos5 and len(FrecRange) m): ss1 = m - - valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1 + + 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 = (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 SA Parameters -------------------------- - + def GetSAParameters(self): #SA en frecuencia pairslist = self.dataOut.groupList num_pairs = len(pairslist) - + vel = self.dataOut.abscissaList - spectra = self.dataOut.data_pre[0] - cspectra = self.dataOut.data_pre[1] - delta_v = vel[1] - vel[0] - + spectra = self.dataOut.data_pre + cspectra = self.dataIn.data_cspc + delta_v = vel[1] - vel[0] + #Calculating the power spectrum spc_pow = numpy.sum(spectra, 3)*delta_v #Normalizing Spectra norm_spectra = spectra/spc_pow #Calculating the norm_spectra at peak - max_spectra = numpy.max(norm_spectra, 3) - + max_spectra = numpy.max(norm_spectra, 3) + #Normalizing Cross Spectra norm_cspectra = numpy.zeros(cspectra.shape) - + for i in range(num_chan): norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:]) - + max_cspectra = numpy.max(norm_cspectra,2) max_cspectra_index = numpy.argmax(norm_cspectra, 2) - + for i in range(num_pairs): cspc_par[i,:,:] = __calculateMoments(norm_cspectra) #------------------- Get Lags ---------------------------------- - + class SALags(Operation): ''' Function GetMoments() @@ -283,19 +1523,19 @@ class SALags(Operation): self.dataOut.data_SNR self.dataOut.groupList self.dataOut.nChannels - + Affected: self.dataOut.data_param - + ''' - def run(self, dataOut): + def run(self, dataOut): data_acf = dataOut.data_pre[0] data_ccf = dataOut.data_pre[1] normFactor_acf = dataOut.normFactor[0] normFactor_ccf = dataOut.normFactor[1] pairs_acf = dataOut.groupList[0] pairs_ccf = dataOut.groupList[1] - + nHeights = dataOut.nHeights absc = dataOut.abscissaList noise = dataOut.noise @@ -306,97 +1546,97 @@ class SALags(Operation): for l in range(len(pairs_acf)): data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:] - + for l in range(len(pairs_ccf)): data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:] - + dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights)) dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc) dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc) return - + # def __getPairsAutoCorr(self, pairsList, nChannels): -# +# # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan -# -# for l in range(len(pairsList)): +# +# for l in range(len(pairsList)): # firstChannel = pairsList[l][0] # secondChannel = pairsList[l][1] -# -# #Obteniendo pares de Autocorrelacion +# +# #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_acf, data_ccf, lagRange): - + lag0 = data_acf.shape[1]/2 #Funcion de Autocorrelacion mean_acf = stats.nanmean(data_acf, axis = 0) - + #Obtencion Indice de TauCross ind_ccf = data_ccf.argmax(axis = 1) #Obtencion Indice de TauAuto ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int') ccf_lag0 = data_ccf[:,lag0,:] - + for i in range(ccf_lag0.shape[0]): ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0) - + #Obtencion de TauCross y TauAuto tau_ccf = lagRange[ind_ccf] tau_acf = lagRange[ind_acf] - + Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0]) - + tau_ccf[Nan1,Nan2] = numpy.nan tau_acf[Nan1,Nan2] = numpy.nan tau = numpy.vstack((tau_ccf,tau_acf)) - + return tau - + def __calculateLag1Phase(self, data, lagTRange): data1 = stats.nanmean(data, axis = 0) lag1 = numpy.where(lagTRange == 0)[0][0] + 1 phase = numpy.angle(data1[lag1,:]) - + return phase - + class SpectralFitting(Operation): ''' Function GetMoments() - + Input: Output: Variables modified: ''' - - def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None): - - + + def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None): + + if path != None: sys.path.append(path) self.dataOut.library = importlib.import_module(file) - + #To be inserted as a parameter groupArray = numpy.array(groupList) -# groupArray = numpy.array([[0,1],[2,3]]) +# groupArray = numpy.array([[0,1],[2,3]]) self.dataOut.groupList = groupArray - + nGroups = groupArray.shape[0] nChannels = self.dataIn.nChannels nHeights=self.dataIn.heightList.size - + #Parameters Array self.dataOut.data_param = None - + #Set constants constants = self.dataOut.library.setConstants(self.dataIn) self.dataOut.constants = constants @@ -405,24 +1645,24 @@ class SpectralFitting(Operation): ippSeconds = self.dataIn.ippSeconds K = self.dataIn.nIncohInt pairsArray = numpy.array(self.dataIn.pairsList) - + #List of possible combinations listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) indCross = numpy.zeros(len(list(listComb)), dtype = 'int') - + if getSNR: listChannels = groupArray.reshape((groupArray.size)) listChannels.sort() noise = self.dataIn.getNoise() self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels]) - - for i in range(nGroups): + + for i in range(nGroups): coord = groupArray[i,:] - + #Input data array data = self.dataIn.data_spc[coord,:,:]/(M*N) data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) - + #Cross Spectra data array for Covariance Matrixes ind = 0 for pairs in listComb: @@ -431,10 +1671,10 @@ class SpectralFitting(Operation): ind += 1 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N) dataCross = dataCross**2/K - + for h in range(nHeights): # print self.dataOut.heightList[h] - + #Input d = data[:,h] @@ -443,7 +1683,7 @@ class SpectralFitting(Operation): ind = 0 for pairs in listComb: #Coordinates in Covariance Matrix - x = pairs[0] + x = pairs[0] y = pairs[1] #Channel Index S12 = dataCross[ind,:,h] @@ -457,15 +1697,15 @@ class SpectralFitting(Operation): LT=L.T dp = numpy.dot(LT,d) - + #Initial values data_spc = self.dataIn.data_spc[coord,:,h] - + if (h>0)and(error1[3]<5): p0 = self.dataOut.data_param[i,:,h-1] else: p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i)) - + try: #Least Squares minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) @@ -478,30 +1718,30 @@ class SpectralFitting(Operation): minp = p0*numpy.nan error0 = numpy.nan error1 = p0*numpy.nan - + #Save - if self.dataOut.data_param is None: + if self.dataOut.data_param == None: self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan - + self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) self.dataOut.data_param[i,:,h] = minp return - + def __residFunction(self, p, dp, LT, constants): fm = self.dataOut.library.modelFunction(p, constants) fmp=numpy.dot(LT,fm) - + return dp-fmp def __getSNR(self, z, noise): - + avg = numpy.average(z, axis=1) SNR = (avg.T-noise)/noise SNR = SNR.T return SNR - + def __chisq(p,chindex,hindex): #similar to Resid but calculates CHI**2 [LT,d,fm]=setupLTdfm(p,chindex,hindex) @@ -509,50 +1749,53 @@ class SpectralFitting(Operation): fmp=numpy.dot(LT,fm) chisq=numpy.dot((dp-fmp).T,(dp-fmp)) return chisq - + 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 __calculateCosDir(self, elev, azim): zen = (90 - elev)*numpy.pi/180 azim = azim*numpy.pi/180 - cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2))) + cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2))) cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2) - + signX = numpy.sign(numpy.cos(azim)) signY = numpy.sign(numpy.sin(azim)) - + cosDirX = numpy.copysign(cosDirX, signX) cosDirY = numpy.copysign(cosDirY, signY) return cosDirX, cosDirY - + 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: @@ -566,37 +1809,37 @@ class WindProfiler(Operation): listPhi = phi.tolist() maxid = listPhi.index(max(listPhi)) minid = listPhi.index(min(listPhi)) - - rango = range(len(phi)) + + 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]): @@ -604,27 +1847,27 @@ class WindProfiler(Operation): # 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): - + def techniqueDBS(self, kwargs): """ 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 """ velRadial0 = kwargs['velRadial'] heiRang = kwargs['heightList'] SNR0 = kwargs['SNR'] - + if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'): theta_x = numpy.array(kwargs['dirCosx']) theta_y = numpy.array(kwargs['dirCosy']) @@ -632,7 +1875,7 @@ class WindProfiler(Operation): elev = numpy.array(kwargs['elevation']) azim = numpy.array(kwargs['azimuth']) theta_x, theta_y = self.__calculateCosDir(elev, azim) - azimuth = kwargs['correctAzimuth'] + azimuth = kwargs['correctAzimuth'] if kwargs.has_key('horizontalOnly'): horizontalOnly = kwargs['horizontalOnly'] else: horizontalOnly = False @@ -647,22 +1890,22 @@ class WindProfiler(Operation): param = param[arrayChannel,:,:] theta_x = theta_x[arrayChannel] theta_y = theta_y[arrayChannel] - - azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) - heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0) + + azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*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, pairs_ccf, azimuth = None): - + nPairs = len(pairs_ccf) posx = numpy.asarray(posx) posy = numpy.asarray(posy) - + #Rotacion Inversa para alinear con el azimuth if azimuth!= None: azimuth = azimuth*math.pi/180 @@ -671,126 +1914,126 @@ class WindProfiler(Operation): else: posx1 = posx posy1 = posy - + #Calculo de Distancias distx = numpy.zeros(nPairs) disty = numpy.zeros(nPairs) dist = numpy.zeros(nPairs) ang = numpy.zeros(nPairs) - + for i in range(nPairs): distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]] - disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]] + disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]] dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2) ang[i] = numpy.arctan2(disty[i],distx[i]) - + return distx, disty, dist, ang - #Calculo de Matrices + #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] nHeights = tau1.shape[1] - vel = numpy.zeros((nPairs,3,nHeights)) + vel = numpy.zeros((nPairs,3,nHeights)) dist1 = numpy.reshape(dist, (dist.size,1)) - + angCos = numpy.cos(ang) angSin = numpy.sin(ang) - - vel0 = dist1*tau1/(2*tau2**2) + + vel0 = dist1*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)): +# +# for l in range(len(pairsList)): # firstChannel = pairsList[l][0] # secondChannel = pairsList[l][1] -# -# #Obteniendo pares de Autocorrelacion +# +# #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): def techniqueSA(self, kwargs): - - """ + + """ 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 """ position_x = kwargs['positionX'] position_y = kwargs['positionY'] azimuth = kwargs['azimuth'] - + if kwargs.has_key('correctFactor'): correctFactor = kwargs['correctFactor'] else: correctFactor = 1 - + groupList = kwargs['groupList'] pairs_ccf = groupList[1] tau = kwargs['tau'] _lambda = kwargs['_lambda'] - + #Cross Correlation pairs obtained # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, 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 + #Metodo Directo distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth) winds = self.__calculateVelHorDir(dist, tau1, tau2, ang) winds = stats.nanmean(winds, axis=0) @@ -806,100 +2049,100 @@ class WindProfiler(Operation): winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) winds = correctFactor*winds return winds - + def __checkTime(self, currentTime, paramInterval, outputInterval): - + dataTime = currentTime + paramInterval deltaTime = dataTime - self.__initime - + if deltaTime >= outputInterval or deltaTime < 0: self.__dataReady = True - return - - def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax, binkm=2): + 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 - ''' -# print arrayMeteor.shape + ''' +# print arrayMeteor.shape #Settings - nInt = (heightMax - heightMin)/binkm + nInt = (heightMax - heightMin)/2 # print nInt nInt = int(nInt) # print nInt - winds = numpy.zeros((2,nInt))*numpy.nan - + winds = numpy.zeros((2,nInt))*numpy.nan + #Filter errors error = numpy.where(arrayMeteor[:,-1] == 0)[0] finalMeteor = arrayMeteor[error,:] - + #Meteor Histogram finalHeights = finalMeteor[:,2] 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 - + 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[:, 6] zen = meteorAux[:, 4]*numpy.pi/180 azim = meteorAux[:, 3]*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 techniqueNSM_SA(self, **kwargs): metArray = kwargs['metArray'] heightList = kwargs['heightList'] timeList = kwargs['timeList'] - + rx_location = kwargs['rx_location'] groupList = kwargs['groupList'] azimuth = kwargs['azimuth'] dfactor = kwargs['dfactor'] k = kwargs['k'] - + azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth) d = dist*dfactor #Phase calculation metArray1 = self.__getPhaseSlope(metArray, heightList, timeList) - + metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities - + velEst = numpy.zeros((heightList.size,2))*numpy.nan azimuth1 = azimuth1*numpy.pi/180 - + for i in range(heightList.size): h = heightList[i] indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0] @@ -912,71 +2155,71 @@ class WindProfiler(Operation): A = numpy.asmatrix(A) A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose() velHor = numpy.dot(A1,velAux) - + velEst[i,:] = numpy.squeeze(velHor) return velEst - + def __getPhaseSlope(self, metArray, heightList, timeList): meteorList = [] #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2 #Putting back together the meteor matrix utctime = metArray[:,0] uniqueTime = numpy.unique(utctime) - + phaseDerThresh = 0.5 ippSeconds = timeList[1] - timeList[0] sec = numpy.where(timeList>1)[0][0] nPairs = metArray.shape[1] - 6 nHeights = len(heightList) - + for t in uniqueTime: metArray1 = metArray[utctime==t,:] # phaseDerThresh = numpy.pi/4 #reducir Phase thresh tmet = metArray1[:,1].astype(int) hmet = metArray1[:,2].astype(int) - + metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1)) metPhase[:,:] = numpy.nan metPhase[:,hmet,tmet] = metArray1[:,6:].T - + #Delete short trails metBool = ~numpy.isnan(metPhase[0,:,:]) heightVect = numpy.sum(metBool, axis = 1) metBool[heightVect phaseDerThresh)) metPhase[phDerAux] = numpy.nan - + #--------------------------METEOR DETECTION ----------------------------------------- indMet = numpy.where(numpy.any(metBool,axis=1))[0] - + for p in numpy.arange(nPairs): phase = metPhase[p,:,:] phDer = metDer[p,:,:] - + for h in indMet: height = heightList[h] phase1 = phase[h,:] #82 phDer1 = phDer[h,:] - + phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap - + indValid = numpy.where(~numpy.isnan(phase1))[0] initMet = indValid[0] endMet = 0 - + for i in range(len(indValid)-1): - + #Time difference inow = indValid[i] inext = indValid[i+1] idiff = inext - inow #Phase difference - phDiff = numpy.abs(phase1[inext] - phase1[inow]) - + phDiff = numpy.abs(phase1[inext] - phase1[inow]) + if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor sizeTrail = inow - initMet + 1 if sizeTrail>3*sec: #Too short meteors @@ -992,28 +2235,28 @@ class WindProfiler(Operation): vel = slope#*height*1000/(k*d) estAux = numpy.array([utctime,p,height, vel, rsq]) meteorList.append(estAux) - initMet = inext + initMet = inext metArray2 = numpy.array(meteorList) - + return metArray2 - + def __calculateAzimuth1(self, rx_location, pairslist, azimuth0): - + azimuth1 = numpy.zeros(len(pairslist)) dist = numpy.zeros(len(pairslist)) - + for i in range(len(rx_location)): ch0 = pairslist[i][0] ch1 = pairslist[i][1] - + diffX = rx_location[ch0][0] - rx_location[ch1][0] diffY = rx_location[ch0][1] - rx_location[ch1][1] azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi dist[i] = numpy.sqrt(diffX**2 + diffY**2) - + azimuth1 -= azimuth0 return azimuth1, dist - + def techniqueNSM_DBS(self, **kwargs): metArray = kwargs['metArray'] heightList = kwargs['heightList'] @@ -1068,39 +2311,39 @@ class WindProfiler(Operation): # noise = dataOut.noise heightList = dataOut.heightList SNR = dataOut.data_SNR - + if technique == 'DBS': - - kwargs['velRadial'] = param[:,1,:] #Radial velocity + + kwargs['velRadial'] = param[:,1,:] #Radial velocity kwargs['heightList'] = heightList kwargs['SNR'] = SNR - + dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.paramInterval - + 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 -# +# pairs = None +# # if kwargs.has_key('correctFactor'): # correctFactor = kwargs['correctFactor'] # else: # correctFactor = 1 - + # tau = dataOut.data_param # _lambda = dataOut.C/dataOut.frequency # pairsList = dataOut.groupList # nChannels = dataOut.nChannels - + kwargs['groupList'] = dataOut.groupList kwargs['tau'] = dataOut.data_param kwargs['_lambda'] = dataOut.C/dataOut.frequency @@ -1108,35 +2351,30 @@ class WindProfiler(Operation): dataOut.data_output = self.techniqueSA(kwargs) dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.timeInterval - - elif technique == 'Meteors': + + elif technique == 'Meteors': dataOut.flagNoData = True self.__dataReady = False - + if kwargs.has_key('nHours'): nHours = kwargs['nHours'] - else: + 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 - - if kwargs.has_key('BinKm'): - binkm = kwargs['BinKm'] - else: - binkm = 2 - + dataOut.outputInterval = nHours*3600 - + if self.__isConfig == False: # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) #Get Initial LTC time @@ -1144,29 +2382,29 @@ class WindProfiler(Operation): self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() self.__isConfig = True - - if self.__buffer is None: + + 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.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready - + if self.__dataReady: dataOut.utctimeInit = self.__initime - + self.__initime += dataOut.outputInterval #to erase time offset - - dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax, binkm) + + dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) dataOut.flagNoData = False self.__buffer = None - + elif technique == 'Meteors1': dataOut.flagNoData = True self.__dataReady = False - + if kwargs.has_key('nMins'): nMins = kwargs['nMins'] else: nMins = 20 @@ -1187,17 +2425,17 @@ class WindProfiler(Operation): else: mode = 'SA' #Borrar luego esto - if dataOut.groupList is None: + if dataOut.groupList == None: dataOut.groupList = [(0,1),(0,2),(1,2)] groupList = dataOut.groupList C = 3e8 freq = 50e6 lamb = C/freq k = 2*numpy.pi/lamb - + timeList = dataOut.abscissaList heightList = dataOut.heightList - + if self.__isConfig == False: dataOut.outputInterval = nMins*60 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) @@ -1208,20 +2446,20 @@ class WindProfiler(Operation): self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() self.__isConfig = True - - if self.__buffer is None: + + 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.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready - + if self.__dataReady: dataOut.utctimeInit = self.__initime self.__initime += dataOut.outputInterval #to erase time offset - + metArray = self.__buffer if mode == 'SA': dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList) @@ -1232,69 +2470,71 @@ class WindProfiler(Operation): self.__buffer = None return - + class EWDriftsEstimation(Operation): - - + + def __init__(self): + Operation.__init__(self) + def __correctValues(self, heiRang, phi, velRadial, SNR): listPhi = phi.tolist() maxid = listPhi.index(max(listPhi)) minid = listPhi.index(min(listPhi)) - - rango = range(len(phi)) + + rango = range(len(phi)) # rango = numpy.delete(rango,maxid) - + heiRang1 = heiRang*math.cos(phi[maxid]) heiRangAux = heiRang*math.cos(phi[minid]) indOut = (heiRang1 < heiRangAux[0]).nonzero() heiRang1 = numpy.delete(heiRang1,indOut) - + velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) SNR1 = numpy.zeros([len(phi),len(heiRang1)]) - + for i in rango: x = heiRang*math.cos(phi[i]) y1 = velRadial[i,:] f1 = interpolate.interp1d(x,y1,kind = 'cubic') - + x1 = heiRang1 y11 = f1(x1) - + y2 = SNR[i,:] f2 = interpolate.interp1d(x,y2,kind = 'cubic') y21 = f2(x1) - + velRadial1[i,:] = y11 SNR1[i,:] = y21 - + return heiRang1, velRadial1, SNR1 def run(self, dataOut, zenith, zenithCorrection): heiRang = dataOut.heightList velRadial = dataOut.data_param[:,3,:] SNR = dataOut.data_SNR - + zenith = numpy.array(zenith) - zenith -= zenithCorrection + zenith -= zenithCorrection zenith *= numpy.pi/180 - + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR) - + alp = zenith[0] bet = zenith[1] - + w_w = velRadial1[0,:] w_e = velRadial1[1,:] - - w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp)) - u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp)) - + + w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp)) + u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp)) + winds = numpy.vstack((u,w)) - + dataOut.heightList = heiRang1 dataOut.data_output = winds dataOut.data_SNR = SNR1 - + dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.timeInterval return @@ -1333,7 +2573,7 @@ class NonSpecularMeteorDetection(Operation): SNR[i] = (power[i]-noise[i])/noise[i] SNRm = numpy.nanmean(SNR, axis = 0) SNRdB = 10*numpy.log10(SNR) - + if mode == 'SA': dataOut.groupList = dataOut.groupList[1] nPairs = data_ccf.shape[0] @@ -1341,22 +2581,22 @@ class NonSpecularMeteorDetection(Operation): phase = numpy.zeros(data_ccf[:,0,:,:].shape) # phase1 = numpy.copy(phase) coh1 = numpy.zeros(data_ccf[:,0,:,:].shape) - + for p in range(nPairs): ch0 = pairsList[p][0] ch1 = pairsList[p][1] ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:]) - phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter -# phase1[p,:,:] = numpy.angle(ccf) #median filter - coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter -# coh1[p,:,:] = numpy.abs(ccf) #median filter + phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter +# phase1[p,:,:] = numpy.angle(ccf) #median filter + coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter +# coh1[p,:,:] = numpy.abs(ccf) #median filter coh = numpy.nanmax(coh1, axis = 0) # struc = numpy.ones((5,1)) # coh = ndimage.morphology.grey_dilation(coh, size=(10,1)) #---------------------- Radial Velocity ---------------------------- phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0) velRad = phaseAux*lamb/(4*numpy.pi*tSamp) - + if allData: boolMetFin = ~numpy.isnan(SNRm) # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0) @@ -1364,31 +2604,31 @@ class NonSpecularMeteorDetection(Operation): #------------------------ Meteor mask --------------------------------- # #SNR mask # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB)) -# +# # #Erase small objects -# boolMet1 = self.__erase_small(boolMet, 2*sec, 5) -# +# boolMet1 = self.__erase_small(boolMet, 2*sec, 5) +# # auxEEJ = numpy.sum(boolMet1,axis=0) # indOver = auxEEJ>nProfiles*0.8 #Use this later # indEEJ = numpy.where(indOver)[0] # indNEEJ = numpy.where(~indOver)[0] -# +# # boolMetFin = boolMet1 -# +# # if indEEJ.size > 0: -# boolMet1[:,indEEJ] = False #Erase heights with EEJ -# +# boolMet1[:,indEEJ] = False #Erase heights with EEJ +# # boolMet2 = coh > cohThresh # boolMet2 = self.__erase_small(boolMet2, 2*sec,5) -# +# # #Final Meteor mask # boolMetFin = boolMet1|boolMet2 - + #Coherence mask boolMet1 = coh > 0.75 struc = numpy.ones((30,1)) boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc) - + #Derivative mask derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0) boolMet2 = derPhase < 0.2 @@ -1405,7 +2645,7 @@ class NonSpecularMeteorDetection(Operation): tmet = coordMet[0] hmet = coordMet[1] - + data_param = numpy.zeros((tmet.size, 6 + nPairs)) data_param[:,0] = utctime data_param[:,1] = tmet @@ -1414,7 +2654,7 @@ class NonSpecularMeteorDetection(Operation): data_param[:,4] = velRad[tmet,hmet] data_param[:,5] = coh[tmet,hmet] data_param[:,6:] = phase[:,tmet,hmet].T - + elif mode == 'DBS': dataOut.groupList = numpy.arange(nChannels) @@ -1422,7 +2662,7 @@ class NonSpecularMeteorDetection(Operation): phase = numpy.angle(data_acf[:,1,:,:]) # phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1)) velRad = phase*lamb/(4*numpy.pi*tSamp) - + #Spectral width # acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1)) # acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1)) @@ -1437,24 +2677,24 @@ class NonSpecularMeteorDetection(Operation): #SNR boolMet1 = (SNRdB>SNRthresh) #SNR mask boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5)) - + #Radial velocity boolMet2 = numpy.abs(velRad) < 20 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5)) - + #Spectral Width boolMet3 = spcWidth < 30 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5)) # boolMetFin = self.__erase_small(boolMet1, 10,5) boolMetFin = boolMet1&boolMet2&boolMet3 - + #Creating data_param coordMet = numpy.where(boolMetFin) cmet = coordMet[0] tmet = coordMet[1] hmet = coordMet[2] - + data_param = numpy.zeros((tmet.size, 7)) data_param[:,0] = utctime data_param[:,1] = cmet @@ -1463,7 +2703,7 @@ class NonSpecularMeteorDetection(Operation): data_param[:,4] = SNR[cmet,tmet,hmet].T data_param[:,5] = velRad[cmet,tmet,hmet].T data_param[:,6] = spcWidth[cmet,tmet,hmet].T - + # self.dataOut.data_param = data_int if len(data_param) == 0: dataOut.flagNoData = True @@ -1473,21 +2713,21 @@ class NonSpecularMeteorDetection(Operation): def __erase_small(self, binArray, threshX, threshY): labarray, numfeat = ndimage.measurements.label(binArray) binArray1 = numpy.copy(binArray) - + for i in range(1,numfeat + 1): auxBin = (labarray==i) auxSize = auxBin.sum() - + x,y = numpy.where(auxBin) widthX = x.max() - x.min() widthY = y.max() - y.min() - + #width X: 3 seg -> 12.5*3 - #width Y: - + #width Y: + if (auxSize < 50) or (widthX < threshX) or (widthY < threshY): binArray1[auxBin] = False - + return binArray1 #--------------- Specular Meteor ---------------- @@ -1497,36 +2737,36 @@ class SMDetection(Operation): 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 - + 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 @@ -1545,9 +2785,9 @@ class SMDetection(Operation): 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): Utc Time | Range Height @@ -1555,21 +2795,21 @@ class SMDetection(Operation): VelRad errorVelRad Phase0 Phase1 Phase2 Phase3 TypeError - - ''' - + + ''' + def run(self, dataOut, hei_ref = None, tauindex = 0, phaseOffsets = None, - cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, + cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, noise_timeStep = 4, noise_multiple = 4, multDet_timeLimit = 1, multDet_rangeLimit = 3, phaseThresh = 20, SNRThresh = 5, hmin = 50, hmax=150, azimuth = 0, channelPositions = None) : - - + + #Getting Pairslist - if channelPositions is None: + if channelPositions == None: # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella meteorOps = SMOperations() @@ -1577,53 +2817,53 @@ class SMDetection(Operation): heiRang = dataOut.getHeiRange() #Get Beacon signal - No Beacon signal anymore # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex]) -# +# # if hei_ref != None: # newheis = numpy.where(self.dataOut.heightList>hei_ref) -# - - +# + + #****************REMOVING HARDWARE PHASE DIFFERENCES*************** # see if the user put in pre defined phase shifts voltsPShift = dataOut.data_pre.copy() - + # if predefinedPhaseShifts != None: # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180 -# +# # # elif beaconPhaseShifts: # # #get hardware phase shifts using beacon signal # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10) # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0) -# +# # else: -# hardwarePhaseShifts = numpy.zeros(5) -# +# hardwarePhaseShifts = numpy.zeros(5) +# # 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 + + #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, dataOut.timeInterval, pairslist0, 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, dataOut.timeInterval) @@ -1633,7 +2873,7 @@ class SMDetection(Operation): #Meteor echoes detection listMeteors = self.__findMeteors(powerNet, signalThresh) #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION ********** - + #************** REMOVE MULTIPLE DETECTIONS (3.5) *************************** #Parameters heiRange = dataOut.getHeiRange() @@ -1643,7 +2883,7 @@ class SMDetection(Operation): #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 @@ -1654,40 +2894,40 @@ class SMDetection(Operation): #Estimation of decay times (Errors N 7, 8, 11) listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, 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, pairslist0, dataOut.timeInterval) + listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval) if len(listMeteors4) > 0: #Setting New Array date = dataOut.utctime arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) - + #Correcting phase offset if phaseOffsets != None: phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets) - + #Second Pairslist pairsList = [] pairx = (0,1) pairy = (2,3) pairsList.append(pairx) pairsList.append(pairy) - + jph = numpy.array([0,0,0,0]) h = (hmin,hmax) arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph) - + # #Calculate AOA (Error N 3, 4) # #JONES ET AL. 1998 # error = arrayParameters[:,-1] # AOAthresh = numpy.pi/8 # phases = -arrayParameters[:,9:13] # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth) -# +# # #Calculate Heights (Error N 13 and 14) # error = arrayParameters[:,-1] # Ranges = arrayParameters[:,2] @@ -1695,73 +2935,73 @@ class SMDetection(Operation): # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax) # error = arrayParameters[:,-1] #********************* END OF PARAMETERS CALCULATION ************************** - - #***************************+ PASS DATA TO NEXT STEP ********************** + + #***************************+ PASS DATA TO NEXT STEP ********************** # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1])) dataOut.data_param = arrayParameters - - if arrayParameters is None: + + if arrayParameters == None: dataOut.flagNoData = True else: dataOut.flagNoData = True - + 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() + 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)) + ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor)) phaseOffset[ind] = numpy.nan - phaseOffset = stats.nanmean(phaseOffset, axis=1) - + 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) + 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: + + 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] @@ -1771,16 +3011,16 @@ class SMDetection(Operation): #Phase Differences phaseDiff = phaseInt - phaseCCF[:,2,:] phaseArrival = phaseInt.reshape(phaseInt.size) - + #Dealias phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival)) # 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 @@ -1790,25 +3030,25 @@ class SMDetection(Operation): 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] - + 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 @@ -1820,21 +3060,21 @@ class SMDetection(Operation): 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] + 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]] - + volts2 = volts[pairslist[i][1]] + for t in range(len(laglist)): - idxT = laglist[t] + idxT = laglist[t] if idxT >= 0: vStacked = numpy.vstack((volts2[idxT:,:], numpy.zeros((idxT, nHeights),dtype='complex'))) @@ -1842,10 +3082,10 @@ class SMDetection(Operation): 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) @@ -1854,133 +3094,133 @@ class SMDetection(Operation): 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] - + 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 - + noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux + noiseAux1 = numpy.mean(arrayBlock) - noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1 - + 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] if isinstance(indUPthresh[j], (int, float)) else indUPthresh[j][0] ##CHECK!!!! - + 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) + + 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,4)) - pairslist1.append((1,3)) + pairslist1.append((0,1)) + pairslist1.append((3,4)) listMeteors1 = [] listPowerSeries = [] listVoltageSeries = [] #volts has the war data - - if frequency == 30.175e6: + + if frequency == 30e6: timeLag = 45*10**-3 else: timeLag = 15*10**-3 - lag = int(numpy.ceil(timeLag/timeInterval)) - + 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 = int(listMeteors[i][0]) - mStart = int(listMeteors[i][1]) - mPeak = int(listMeteors[i][2]) - mEnd = int(listMeteors[i][3]) - + 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] @@ -1988,15 +3228,15 @@ class SMDetection(Operation): 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 = 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: @@ -2004,7 +3244,7 @@ class SMDetection(Operation): 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] @@ -2013,17 +3253,17 @@ class SMDetection(Operation): 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: + 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) @@ -2038,50 +3278,50 @@ class SMDetection(Operation): #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: + + + 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 - + + return listMeteors1, listPowerSeries, listVoltageSeries + def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency): - + threshError = 10 #Depending if it is 30 or 50 MHz - if frequency == 30.175e6: + if frequency == 30e6: timeLag = 45*10**-3 else: timeLag = 15*10**-3 - lag = int(numpy.ceil(timeLag/timeInterval)) - + lag = numpy.ceil(timeLag/timeInterval) + listMeteors1 = [] - + for i in range(len(listMeteors)): meteorPower = listPower[i] meteorAux = listMeteors[i] - + if meteorAux[-1] == 0: - try: + try: indmax = meteorPower.argmax() indlag = indmax + lag - + y = meteorPower[indlag:] x = numpy.arange(0, y.size)*timeLag - + #first guess a = y[0] tau = timeLag @@ -2090,26 +3330,26 @@ class SMDetection(Operation): 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 + 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 - - + meteorAux[-1] = 11 + + except: - meteorAux[-1] = 11 - - + meteorAux[-1] = 11 + + listMeteors1.append(meteorAux) - + return listMeteors1 #Exponential Function @@ -2117,45 +3357,45 @@ class SMDetection(Operation): 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,4)) - pairslist1.append((1,3)) + 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 = 30.175e6 - + freq = 30e6 + listMeteors1 = [] - + for i in range(len(listMeteors)): meteorAux = listMeteors[i] if meteorAux[-1] == 0: mStart = listMeteors[i][1] - mPeak = listMeteors[i][2] + 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,4,2,3],0]) - + 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: + + 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 @@ -2164,51 +3404,51 @@ class SMDetection(Operation): 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: + 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), 13)) - + #Date inclusion # date = re.findall(r'\((.*?)\)', date) # date = date[0].split(',') # date = map(int, date) -# +# # if len(date)<6: # date.append(0) -# +# # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]] # arrayDate = numpy.tile(date, (len(listMeteors), 1)) arrayDate = numpy.tile(date, (len(listMeteors))) - + #Meteor array # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) - + #Parameters Array arrayParameters[:,0] = arrayDate #Date arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range @@ -2216,13 +3456,13 @@ class SMDetection(Operation): arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases arrayParameters[:,-1] = arrayMeteors[:,-1] #Error - + return arrayParameters - + class CorrectSMPhases(Operation): - + def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None): - + arrayParameters = dataOut.data_param pairsList = [] pairx = (0,1) @@ -2230,49 +3470,49 @@ class CorrectSMPhases(Operation): pairsList.append(pairx) pairsList.append(pairy) jph = numpy.zeros(4) - + phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets) arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets))) - + meteorOps = SMOperations() - if channelPositions is None: + if channelPositions == None: # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella - + pairslist0, distances = meteorOps.getPhasePairs(channelPositions) h = (hmin,hmax) - + arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph) - + dataOut.data_param = arrayParameters return class SMPhaseCalibration(Operation): - + __buffer = None __initime = None __dataReady = False - + __isConfig = False - + def __checkTime(self, currentTime, initTime, paramInterval, outputInterval): - + dataTime = currentTime + paramInterval deltaTime = dataTime - initTime - + if deltaTime >= outputInterval or deltaTime < 0: return True - + return False - + def __getGammas(self, pairs, d, phases): gammas = numpy.zeros(2) - + for i in range(len(pairs)): - + pairi = pairs[i] phip3 = phases[:,pairi[0]] @@ -2286,7 +3526,7 @@ class SMPhaseCalibration(Operation): jgamma = numpy.angle(numpy.exp(1j*jgamma)) # jgamma[jgamma>numpy.pi] -= 2*numpy.pi # jgamma[jgamma<-numpy.pi] += 2*numpy.pi - + #Revised distribution jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi)) @@ -2295,39 +3535,39 @@ class SMPhaseCalibration(Operation): rmin = -0.5*numpy.pi rmax = 0.5*numpy.pi phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax)) - + meteorsY = phaseHisto[0] phasesX = phaseHisto[1][:-1] width = phasesX[1] - phasesX[0] phasesX += width/2 - + #Gaussian aproximation bpeak = meteorsY.argmax() peak = meteorsY.max() jmin = bpeak - 5 jmax = bpeak + 5 + 1 - + if jmin<0: jmin = 0 jmax = 6 elif jmax > meteorsY.size: jmin = meteorsY.size - 6 jmax = meteorsY.size - + x0 = numpy.array([peak,bpeak,50]) coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax])) - + #Gammas gammas[i] = coeff[0][1] - + return gammas - + def __residualFunction(self, coeffs, y, t): - + return y - self.__gauss_function(t, coeffs) def __gauss_function(self, t, coeffs): - + return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2) def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray): @@ -2348,16 +3588,16 @@ class SMPhaseCalibration(Operation): max_xangle = range_angle[iz]/2 + center_xangle min_yangle = -range_angle[iz]/2 + center_yangle max_yangle = range_angle[iz]/2 + center_yangle - + inc_x = (max_xangle-min_xangle)/nstepsx inc_y = (max_yangle-min_yangle)/nstepsy - + alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle penalty = numpy.zeros((nstepsx,nstepsy)) jph_array = numpy.zeros((nchan,nstepsx,nstepsy)) jph = numpy.zeros(nchan) - + # Iterations looking for the offset for iy in range(int(nstepsy)): for ix in range(int(nstepsx)): @@ -2387,24 +3627,24 @@ class SMPhaseCalibration(Operation): error = meteorsArray1[:,-1] ind1 = numpy.where(error==0)[0] penalty[ix,iy] = ind1.size - + i,j = numpy.unravel_index(penalty.argmax(), penalty.shape) phOffset = jph_array[:,i,j] - + center_xangle = phOffset[pairx[1]] center_yangle = phOffset[pairy[1]] - + phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j])) - phOffset = phOffset*180/numpy.pi + phOffset = phOffset*180/numpy.pi return phOffset - - + + def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1): - + dataOut.flagNoData = True - self.__dataReady = False + self.__dataReady = False dataOut.outputInterval = nHours*3600 - + if self.__isConfig == False: # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) #Get Initial LTC time @@ -2412,19 +3652,19 @@ class SMPhaseCalibration(Operation): self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() self.__isConfig = True - - if self.__buffer is None: + + if self.__buffer == None: self.__buffer = dataOut.data_param.copy() else: self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) - + self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready - + if self.__dataReady: dataOut.utctimeInit = self.__initime self.__initime += dataOut.outputInterval #to erase time offset - + freq = dataOut.frequency c = dataOut.C #m/s lamb = c/freq @@ -2452,7 +3692,7 @@ class SMPhaseCalibration(Operation): else: pairs.append((2,3)) # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb] - + meteorsArray = self.__buffer error = meteorsArray[:,-1] boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14) @@ -2460,7 +3700,7 @@ class SMPhaseCalibration(Operation): meteorsArray = meteorsArray[ind1,:] meteorsArray[:,-1] = 0 phases = meteorsArray[:,8:12] - + #Calculate Gammas gammas = self.__getGammas(pairs, distances, phases) # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180 @@ -2469,24 +3709,23 @@ class SMPhaseCalibration(Operation): phasesOff = phasesOff.reshape((1,phasesOff.size)) dataOut.data_output = -phasesOff dataOut.flagNoData = False - dataOut.channelList = pairslist0 self.__buffer = None - - + + return - + class SMOperations(): - + def __init__(self): - + return - + def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph): - + arrayParameters = arrayParameters0.copy() hmin = h[0] hmax = h[1] - + #Calculate AOA (Error N 3, 4) #JONES ET AL. 1998 AOAthresh = numpy.pi/8 @@ -2494,72 +3733,72 @@ class SMOperations(): phases = -arrayParameters[:,8:12] + jph # phases = numpy.unwrap(phases) arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth) - + #Calculate Heights (Error N 13 and 14) error = arrayParameters[:,-1] Ranges = arrayParameters[:,1] zenith = arrayParameters[:,4] arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax) - + #----------------------- Get Final data ------------------------------------ # error = arrayParameters[:,-1] # ind1 = numpy.where(error==0)[0] # arrayParameters = arrayParameters[ind1,:] - + return arrayParameters - + def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth): - + arrayAOA = numpy.zeros((phases.shape[0],3)) cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions) - + 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 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0] error[indError] = 0 #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 + 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 + error[indInvalid] = 4 return arrayAOA, error - + def __getDirectionCosines(self, arrayPhase, pairsList, distances): - + #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): ph0 = arrayPhase[:,pairsList[i][0]] ph1 = arrayPhase[:,pairsList[i][1]] d0 = distances[pairsList[i][0]] d1 = distances[pairsList[i][1]] - - ph0_aux = ph0 + ph1 + + ph0_aux = ph0 + ph1 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux)) # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi -# ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi +# ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi #First Estimation cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1)) - + #Most-Accurate Second Estimation phi1_aux = ph0 - ph1 phi1_aux = phi1_aux.reshape(phi1_aux.size,1) #Direction Cosine 1 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1)) - + #Searching the correct Direction Cosine cosdir0_aux = cosdir0[:,i] cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1) @@ -2568,59 +3807,59 @@ class SMOperations(): 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][:len(ind_h)] 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 + #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0] error[indError] = 0 - indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0] + 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 - + error[indInvalid1] = 13 + return heights, error - + def getPhasePairs(self, channelPositions): chanPos = numpy.array(channelPositions) listOper = list(itertools.combinations(range(5),2)) - + distances = numpy.zeros(4) axisX = [] axisY = [] @@ -2628,15 +3867,15 @@ class SMOperations(): distY = numpy.zeros(3) ix = 0 iy = 0 - + pairX = numpy.zeros((2,2)) pairY = numpy.zeros((2,2)) - + for i in range(len(listOper)): pairi = listOper[i] - + posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:]) - + if posDif[0] == 0: axisY.append(pairi) distY[iy] = posDif[1] @@ -2645,7 +3884,7 @@ class SMOperations(): axisX.append(pairi) distX[ix] = posDif[0] ix += 1 - + for i in range(2): if i==0: dist0 = distX @@ -2653,7 +3892,7 @@ class SMOperations(): else: dist0 = distY axis0 = axisY - + side = numpy.argsort(dist0)[:-1] axis0 = numpy.array(axis0)[side,:] chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0]) @@ -2661,7 +3900,7 @@ class SMOperations(): side = axis1[axis1 != chanC] diff1 = chanPos[chanC,i] - chanPos[side[0],i] diff2 = chanPos[chanC,i] - chanPos[side[1],i] - if diff1<0: + if diff1<0: chan2 = side[0] d2 = numpy.abs(diff1) chan1 = side[1] @@ -2671,7 +3910,7 @@ class SMOperations(): d2 = numpy.abs(diff2) chan1 = side[0] d1 = numpy.abs(diff1) - + if i==0: chanCX = chanC chan1X = chan1 @@ -2683,10 +3922,10 @@ class SMOperations(): chan2Y = chan2 distances[2:4] = numpy.array([d1,d2]) # axisXsides = numpy.reshape(axisX[ix,:],4) -# +# # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0]) # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0]) -# +# # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0] # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0] # channel25X = int(pairX[0,ind25X]) @@ -2695,59 +3934,59 @@ class SMOperations(): # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0] # channel25Y = int(pairY[0,ind25Y]) # channel20Y = int(pairY[1,ind20Y]) - + # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)] - pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)] - + pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)] + return pairslist, distances # 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 +# 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 +# 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 +# phi0_aux[indcsi] -= 2*numpy.pi # indcsi = numpy.where(phi0_aux < -numpy.pi) -# phi0_aux[indcsi] += 2*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) @@ -2756,50 +3995,51 @@ class SMOperations(): # 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] +# +# 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] +# #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 +# error[indInvalid1] = 13 +# +# return heights, error + \ No newline at end of file diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 0ba384e..299b3d4 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -1,3 +1,5 @@ +import itertools + import numpy from jroproc_base import ProcessingUnit, Operation @@ -109,7 +111,10 @@ class SpectraProc(ProcessingUnit): if self.dataIn.type == "Spectra": self.dataOut.copy(self.dataIn) -# self.__selectPairs(pairsList) + if not pairsList: + pairsList = itertools.combinations(self.dataOut.channelList, 2) + if self.dataOut.data_cspc is not None: + self.__selectPairs(pairsList) return True if self.dataIn.type == "Voltage": @@ -178,27 +183,21 @@ class SpectraProc(ProcessingUnit): def __selectPairs(self, pairsList): - if channelList == None: + if not pairsList: return - pairsIndexListSelected = [] - - for thisPair in pairsList: + pairs = [] + pairsIndex = [] - if thisPair not in self.dataOut.pairsList: + for pair in pairsList: + if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList: continue - - pairIndex = self.dataOut.pairsList.index(thisPair) - - pairsIndexListSelected.append(pairIndex) - - if not pairsIndexListSelected: - self.dataOut.data_cspc = None - self.dataOut.pairsList = [] - return - - self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected] - self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected] + pairs.append(pair) + pairsIndex.append(pairs.index(pair)) + + self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex] + self.dataOut.pairsList = pairs + self.dataOut.pairsIndexList = pairsIndex return diff --git a/schainpy/model/utils/jroutils_publish.py b/schainpy/model/utils/jroutils_publish.py index b4ad218..726065a 100644 --- a/schainpy/model/utils/jroutils_publish.py +++ b/schainpy/model/utils/jroutils_publish.py @@ -15,6 +15,7 @@ from multiprocessing import Process from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit from schainpy.model.data.jrodata import JROData +from schainpy.utils import log MAXNUMX = 100 MAXNUMY = 100 @@ -30,14 +31,13 @@ def roundFloats(obj): return round(obj, 2) def decimate(z, MAXNUMY): - # dx = int(len(self.x)/self.__MAXNUMX) + 1 - dy = int(len(z[0])/MAXNUMY) + 1 return z[::, ::dy] class throttle(object): - """Decorator that prevents a function from being called more than once every + ''' + Decorator that prevents a function from being called more than once every time period. To create a function that cannot be called more than once a minute, but will sleep until it can be called: @@ -48,7 +48,7 @@ class throttle(object): for i in range(10): foo() print "This function has run %s times." % i - """ + ''' def __init__(self, seconds=0, minutes=0, hours=0): self.throttle_period = datetime.timedelta( @@ -72,9 +72,169 @@ class throttle(object): return wrapper +class Data(object): + ''' + Object to hold data to be plotted + ''' + + def __init__(self, plottypes, throttle_value): + self.plottypes = plottypes + self.throttle = throttle_value + self.ended = False + self.__times = [] + + def __str__(self): + dum = ['{}{}'.format(key, self.shape(key)) for key in self.data] + return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times)) + + def __len__(self): + return len(self.__times) + + def __getitem__(self, key): + if key not in self.data: + raise KeyError(log.error('Missing key: {}'.format(key))) + + if 'spc' in key: + ret = self.data[key] + else: + ret = numpy.array([self.data[key][x] for x in self.times]) + if ret.ndim > 1: + ret = numpy.swapaxes(ret, 0, 1) + return ret + + def setup(self): + ''' + Configure object + ''' + + self.ended = False + self.data = {} + self.__times = [] + self.__heights = [] + self.__all_heights = set() + for plot in self.plottypes: + self.data[plot] = {} + + def shape(self, key): + ''' + Get the shape of the one-element data for the given key + ''' + + if len(self.data[key]): + if 'spc' in key: + return self.data[key].shape + return self.data[key][self.__times[0]].shape + return (0,) + + def update(self, dataOut): + ''' + Update data object with new dataOut + ''' + + tm = dataOut.utctime + if tm in self.__times: + return + + self.parameters = getattr(dataOut, 'parameters', []) + self.pairs = dataOut.pairsList + self.channels = dataOut.channelList + self.xrange = (dataOut.getFreqRange(1)/1000. , dataOut.getAcfRange(1) , dataOut.getVelRange(1)) + self.interval = dataOut.getTimeInterval() + self.__heights.append(dataOut.heightList) + self.__all_heights.update(dataOut.heightList) + self.__times.append(tm) + + for plot in self.plottypes: + if plot == 'spc': + z = dataOut.data_spc/dataOut.normFactor + self.data[plot] = 10*numpy.log10(z) + if plot == 'cspc': + self.data[plot] = dataOut.data_cspc + if plot == 'noise': + self.data[plot][tm] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + if plot == 'rti': + self.data[plot][tm] = dataOut.getPower() + if plot == 'snr_db': + self.data['snr'][tm] = dataOut.data_SNR + if plot == 'snr': + self.data[plot][tm] = 10*numpy.log10(dataOut.data_SNR) + if plot == 'dop': + self.data[plot][tm] = 10*numpy.log10(dataOut.data_DOP) + if plot == 'mean': + self.data[plot][tm] = dataOut.data_MEAN + if plot == 'std': + self.data[plot][tm] = dataOut.data_STD + if plot == 'coh': + self.data[plot][tm] = dataOut.getCoherence() + if plot == 'phase': + self.data[plot][tm] = dataOut.getCoherence(phase=True) + if plot == 'output': + self.data[plot][tm] = dataOut.data_output + if plot == 'param': + self.data[plot][tm] = dataOut.data_param + + def normalize_heights(self): + ''' + Ensure same-dimension of the data for different heighList + ''' + + H = numpy.array(list(self.__all_heights)) + H.sort() + for key in self.data: + shape = self.shape(key)[:-1] + H.shape + for tm, obj in self.data[key].items(): + h = self.__heights[self.__times.index(tm)] + if H.size == h.size: + continue + index = numpy.where(numpy.in1d(H, h))[0] + dummy = numpy.zeros(shape) + numpy.nan + if len(shape) == 2: + dummy[:, index] = obj + else: + dummy[index] = obj + self.data[key][tm] = dummy + + self.__heights = [H for tm in self.__times] + + def jsonify(self, decimate=False): + ''' + Convert data to json + ''' + + ret = {} + tm = self.times[-1] + + for key, value in self.data: + if key in ('spc', 'cspc'): + ret[key] = roundFloats(self.data[key].to_list()) + else: + ret[key] = roundFloats(self.data[key][tm].to_list()) + + ret['timestamp'] = tm + ret['interval'] = self.interval + + @property + def times(self): + ''' + Return the list of times of the current data + ''' + + ret = numpy.array(self.__times) + ret.sort() + return ret + + @property + def heights(self): + ''' + Return the list of heights of the current data + ''' + + return numpy.array(self.__heights[-1]) class PublishData(Operation): - """Clase publish.""" + ''' + Operation to send data over zmq. + ''' def __init__(self, **kwargs): """Inicio.""" @@ -86,11 +246,11 @@ class PublishData(Operation): def on_disconnect(self, client, userdata, rc): if rc != 0: - print("Unexpected disconnection.") + log.warning('Unexpected disconnection.') self.connect() def connect(self): - print 'trying to connect' + log.warning('trying to connect') try: self.client.connect( host=self.host, @@ -104,7 +264,7 @@ class PublishData(Operation): # retain=True # ) except: - print "MQTT Conection error." + log.error('MQTT Conection error.') self.client = False def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, verbose=True, **kwargs): @@ -119,8 +279,7 @@ class PublishData(Operation): self.zeromq = zeromq self.mqtt = kwargs.get('plottype', 0) self.client = None - self.verbose = verbose - self.dataOut.firstdata = True + self.verbose = verbose setup = [] if mqtt is 1: self.client = mqtt.Client( @@ -175,7 +334,6 @@ class PublishData(Operation): 'type': self.plottype, 'yData': yData } - # print payload elif self.plottype in ('rti', 'power'): data = getattr(self.dataOut, 'data_spc') @@ -229,15 +387,16 @@ class PublishData(Operation): 'timestamp': 'None', 'type': None } - # print 'Publishing data to {}'.format(self.host) + self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0) if self.zeromq is 1: if self.verbose: - print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime) + log.log( + '{} - {}'.format(self.dataOut.type, self.dataOut.datatime), + 'Sending' + ) self.zmq_socket.send_pyobj(self.dataOut) - self.dataOut.firstdata = False - def run(self, dataOut, **kwargs): self.dataOut = dataOut @@ -252,6 +411,7 @@ class PublishData(Operation): if self.zeromq is 1: self.dataOut.finished = True self.zmq_socket.send_pyobj(self.dataOut) + time.sleep(0.1) self.zmq_socket.close() if self.client: self.client.loop_stop() @@ -280,7 +440,7 @@ class ReceiverData(ProcessingUnit): self.receiver = self.context.socket(zmq.PULL) self.receiver.bind(self.address) time.sleep(0.5) - print '[Starting] ReceiverData from {}'.format(self.address) + log.success('ReceiverData from {}'.format(self.address)) def run(self): @@ -290,8 +450,9 @@ class ReceiverData(ProcessingUnit): self.isConfig = True self.dataOut = self.receiver.recv_pyobj() - print '[Receiving] {} - {}'.format(self.dataOut.type, - self.dataOut.datatime.ctime()) + log.log('{} - {}'.format(self.dataOut.type, + self.dataOut.datatime.ctime(),), + 'Receiving') class PlotterReceiver(ProcessingUnit, Process): @@ -305,7 +466,6 @@ class PlotterReceiver(ProcessingUnit, Process): self.mp = False self.isConfig = False self.isWebConfig = False - self.plottypes = [] self.connections = 0 server = kwargs.get('server', 'zmq.pipe') plot_server = kwargs.get('plot_server', 'zmq.web') @@ -325,19 +485,13 @@ class PlotterReceiver(ProcessingUnit, Process): self.realtime = kwargs.get('realtime', False) self.throttle_value = kwargs.get('throttle', 5) self.sendData = self.initThrottle(self.throttle_value) + self.dates = [] self.setup() def setup(self): - self.data = {} - self.data['times'] = [] - for plottype in self.plottypes: - self.data[plottype] = {} - self.data['noise'] = {} - self.data['throttle'] = self.throttle_value - self.data['ENDED'] = False - self.isConfig = True - self.data_web = {} + self.data = Data(self.plottypes, self.throttle_value) + self.isConfig = True def event_monitor(self, monitor): @@ -354,15 +508,13 @@ class PlotterReceiver(ProcessingUnit, Process): self.connections += 1 if evt['event'] == 512: pass - if self.connections == 0 and self.started is True: - self.ended = True evt.update({'description': events[evt['event']]}) if evt['event'] == zmq.EVENT_MONITOR_STOPPED: break monitor.close() - print("event monitor thread done!") + print('event monitor thread done!') def initThrottle(self, throttle_value): @@ -372,65 +524,16 @@ class PlotterReceiver(ProcessingUnit, Process): return sendDataThrottled - def send(self, data): - # print '[sending] data=%s size=%s' % (data.keys(), len(data['times'])) + log.success('Sending {}'.format(data), self.name) self.sender.send_pyobj(data) - - def update(self): - t = self.dataOut.utctime - - if t in self.data['times']: - return - - self.data['times'].append(t) - self.data['dataOut'] = self.dataOut - - for plottype in self.plottypes: - if plottype == 'spc': - z = self.dataOut.data_spc/self.dataOut.normFactor - self.data[plottype] = 10*numpy.log10(z) - self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor) - if plottype == 'cspc': - jcoherence = self.dataOut.data_cspc/numpy.sqrt(self.dataOut.data_spc*self.dataOut.data_spc) - self.data['cspc_coh'] = numpy.abs(jcoherence) - self.data['cspc_phase'] = numpy.arctan2(jcoherence.imag, jcoherence.real)*180/numpy.pi - if plottype == 'rti': - self.data[plottype][t] = self.dataOut.getPower() - if plottype == 'snr': - self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR) - if plottype == 'dop': - self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP) - if plottype == 'mean': - self.data[plottype][t] = self.dataOut.data_MEAN - if plottype == 'std': - self.data[plottype][t] = self.dataOut.data_STD - if plottype == 'coh': - self.data[plottype][t] = self.dataOut.getCoherence() - if plottype == 'phase': - self.data[plottype][t] = self.dataOut.getCoherence(phase=True) - if plottype == 'output': - self.data[plottype][t] = self.dataOut.data_output - if plottype == 'param': - self.data[plottype][t] = self.dataOut.data_param - if self.realtime: - self.data_web['timestamp'] = t - if plottype == 'spc': - self.data_web[plottype] = roundFloats(decimate(self.data[plottype]).tolist()) - elif plottype == 'cspc': - self.data_web['cspc_coh'] = roundFloats(decimate(self.data['cspc_coh']).tolist()) - self.data_web['cspc_phase'] = roundFloats(decimate(self.data['cspc_phase']).tolist()) - elif plottype == 'noise': - self.data_web['noise'] = roundFloats(self.data['noise'][t].tolist()) - else: - self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist()) - self.data_web['interval'] = self.dataOut.getTimeInterval() - self.data_web['type'] = plottype - def run(self): - print '[Starting] {} from {}'.format(self.name, self.address) + log.success( + 'Starting from {}'.format(self.address), + self.name + ) self.context = zmq.Context() self.receiver = self.context.socket(zmq.PULL) @@ -447,39 +550,39 @@ class PlotterReceiver(ProcessingUnit, Process): else: self.sender.bind("ipc:///tmp/zmq.plots") - time.sleep(3) + time.sleep(2) t = Thread(target=self.event_monitor, args=(monitor,)) t.start() while True: - self.dataOut = self.receiver.recv_pyobj() - # print '[Receiving] {} - {}'.format(self.dataOut.type, - # self.dataOut.datatime.ctime()) - - self.update() + dataOut = self.receiver.recv_pyobj() + dt = datetime.datetime.fromtimestamp(dataOut.utctime).date() + sended = False + if dt not in self.dates: + if self.data: + self.data.ended = True + self.send(self.data) + sended = True + self.data.setup() + self.dates.append(dt) - if self.dataOut.firstdata is True: - self.data['STARTED'] = True + self.data.update(dataOut) - if self.dataOut.finished is True: - self.send(self.data) + if dataOut.finished is True: self.connections -= 1 - if self.connections == 0 and self.started: - self.ended = True - self.data['ENDED'] = True + if self.connections == 0 and dt in self.dates: + self.data.ended = True self.send(self.data) - self.setup() - self.started = False + self.data.setup() else: if self.realtime: self.send(self.data) - self.sender_web.send_string(json.dumps(self.data_web)) + # self.sender_web.send_string(self.data.jsonify()) else: - self.sendData(self.send, self.data) - self.started = True + if not sended: + self.sendData(self.send, self.data) - self.data['STARTED'] = False return def sendToWeb(self): @@ -496,6 +599,6 @@ class PlotterReceiver(ProcessingUnit, Process): time.sleep(1) for kwargs in self.operationKwargs.values(): if 'plot' in kwargs: - print '[Sending] Config data to web for {}'.format(kwargs['code'].upper()) + log.success('[Sending] Config data to web for {}'.format(kwargs['code'].upper())) sender_web_config.send_string(json.dumps(kwargs)) - self.isWebConfig = True + self.isWebConfig = True \ No newline at end of file