schain Commit - r1391:fba03565a781 · Jicamarca Repository

clean Rayleigh funcionando, tiempo de ejecución elevado usando todos los pares

joabAM -

r1391:fba03565a781

parent child

Context file:

r1391:fba03565a781

Collapse all files

The requested changes are too big and content was truncated. Show full diff

schainpy/model/graphics/jroplot_base.py +2 0

             # Copyright (c) 2012-2020 Jicamarca Radio Observatory
             # All rights reserved.
             #
             # Distributed under the terms of the BSD 3-clause license.
             """Base class to create plot operations
             """
             import os
             import sys
             import zmq
             import time
             import numpy
             import datetime
             from collections import deque
             from functools import wraps
             from threading import Thread
             import matplotlib
             if 'BACKEND' in os.environ:
                 matplotlib.use(os.environ['BACKEND'])
             elif 'linux' in sys.platform:
                 matplotlib.use("TkAgg")
             elif 'darwin' in sys.platform:
                 matplotlib.use('MacOSX')
             else:
                 from schainpy.utils import log
                 log.warning('Using default Backend="Agg"', 'INFO')
                 matplotlib.use('Agg')
             import matplotlib.pyplot as plt
             from matplotlib.patches import Polygon
             from mpl_toolkits.axes_grid1 import make_axes_locatable
             from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
             from schainpy.model.data.jrodata import PlotterData
             from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
             from schainpy.utils import log
             jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
             blu_values = matplotlib.pyplot.get_cmap(
                 'seismic_r', 20)(numpy.arange(20))[10:15]
             ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
                 'jro', numpy.vstack((blu_values, jet_values)))
             matplotlib.pyplot.register_cmap(cmap=ncmap)
             CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
                                                'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
             EARTH_RADIUS = 6.3710e3
             def ll2xy(lat1, lon1, lat2, lon2):
                 p = 0.017453292519943295
                 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
                     numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
                 r = 12742 * numpy.arcsin(numpy.sqrt(a))
                 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
                                       * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
                 theta = -theta + numpy.pi/2
                 return r*numpy.cos(theta), r*numpy.sin(theta)
             def km2deg(km):
                 '''
                 Convert distance in km to degrees
                 '''
                 return numpy.rad2deg(km/EARTH_RADIUS)
             def figpause(interval):
                 backend = plt.rcParams['backend']
                 if backend in matplotlib.rcsetup.interactive_bk:
                     figManager = matplotlib._pylab_helpers.Gcf.get_active()
                     if figManager is not None:
                         canvas = figManager.canvas
                         if canvas.figure.stale:
                             canvas.draw()
                         try:
                             canvas.start_event_loop(interval)
                         except:
                             pass
                         return
             def popup(message):
                 '''
                 '''
                 fig = plt.figure(figsize=(12, 8), facecolor='r')
                 text = '\n'.join([s.strip() for s in message.split(':')])
                 fig.text(0.01, 0.5, text, ha='left', va='center',
                          size='20', weight='heavy', color='w')
                 fig.show()
                 figpause(1000)
             class Throttle(object):
                 '''
                 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:
                 @Throttle(minutes=1)
                 def foo():
                   pass
                 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(
                         seconds=seconds, minutes=minutes, hours=hours
                     )
                     self.time_of_last_call = datetime.datetime.min
                 def __call__(self, fn):
                     @wraps(fn)
                     def wrapper(*args, **kwargs):
                         coerce = kwargs.pop('coerce', None)
                         if coerce:
                             self.time_of_last_call = datetime.datetime.now()
                             return fn(*args, **kwargs)
                         else:
                             now = datetime.datetime.now()
                             time_since_last_call = now - self.time_of_last_call
                             time_left = self.throttle_period - time_since_last_call
                             if time_left > datetime.timedelta(seconds=0):
                                 return
                         self.time_of_last_call = datetime.datetime.now()
                         return fn(*args, **kwargs)
                     return wrapper
             def apply_throttle(value):
                 @Throttle(seconds=value)
                 def fnThrottled(fn):
                     fn()
                 return fnThrottled
             @MPDecorator
             class Plot(Operation):
                 """Base class for Schain plotting operations
                 This class should never be use directtly you must subclass a new operation,
                 children classes must be defined as follow:
                 ExamplePlot(Plot):
                     CODE = 'code'
                     colormap = 'jet'
                     plot_type = 'pcolor' # options are ('pcolor', 'pcolorbuffer', 'scatter', 'scatterbuffer')
                     def setup(self):
                         pass
                     def plot(self):
                         pass
                 """
                 CODE = 'Figure'
                 colormap = 'jet'
                 bgcolor = 'white'
                 buffering = True
                 __missing = 1E30
                 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
                              'showprofile']
                 def __init__(self):
                     Operation.__init__(self)
                     self.isConfig = False
                     self.isPlotConfig = False
                     self.save_time = 0
                     self.sender_time = 0
                     self.data = None
                     self.firsttime = True
                     self.sender_queue = deque(maxlen=10)
                     self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
                 def __fmtTime(self, x, pos):
                     '''
                     '''
                     return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
                 def __setup(self, **kwargs):
                     '''
                     Initialize variables
                     '''
                     self.figures = []
                     self.axes = []
                     self.cb_axes = []
+                    self.pf_axes = []
                     self.localtime = kwargs.pop('localtime', True)
                     self.show = kwargs.get('show', True)
                     self.save = kwargs.get('save', False)
                     self.save_period = kwargs.get('save_period', 0)
                     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.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.labels = kwargs.get('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)
                     self.xmax = kwargs.get('xmax', None)
                     self.xrange = kwargs.get('xrange', 12)
                     self.xscale = kwargs.get('xscale', None)
                     self.ymin = kwargs.get('ymin', None)
                     self.ymax = kwargs.get('ymax', None)
                     self.yscale = kwargs.get('yscale', None)
                     self.xlabel = kwargs.get('xlabel', None)
                     self.attr_time = kwargs.get('attr_time', 'utctime')
                     self.attr_data = kwargs.get('attr_data', 'data_param')
                     self.decimation = kwargs.get('decimation', None)
                     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.channels = kwargs.get('channels', None)
                     self.titles = kwargs.get('titles', [])
                     self.polar = False
                     self.type = kwargs.get('type', 'iq')
                     self.grid = kwargs.get('grid', False)
                     self.pause = kwargs.get('pause', False)
                     self.save_code = kwargs.get('save_code', self.CODE)
                     self.throttle = kwargs.get('throttle', 0)
                     self.exp_code = kwargs.get('exp_code', None)
                     self.server = kwargs.get('server', False)
                     self.sender_period = kwargs.get('sender_period', 60)
                     self.tag = kwargs.get('tag', '')
                     self.height_index = kwargs.get('height_index', None)
                     self.__throttle_plot = apply_throttle(self.throttle)
                     code = self.attr_data if self.attr_data else self.CODE
                     self.data = PlotterData(self.CODE, self.exp_code, self.localtime)
+                    self.tmin = kwargs.get('tmin', None)
                     if self.server:
                         if not self.server.startswith('tcp://'):
                             self.server = 'tcp://{}'.format(self.server)
                         log.success(
                             'Sending to server: {}'.format(self.server),
                             self.name
                         )
                     if isinstance(self.attr_data, str):
                         self.attr_data = [self.attr_data]
                 def __setup_plot(self):
                     '''
                     Common setup for all figures, here figures and axes are created
                     '''
                     self.setup()
                     self.time_label = 'LT' if self.localtime else 'UTC'
                     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, polar=self.polar)
                             ax.tick_params(labelsize=8)
                             ax.firsttime = True
                             ax.index = 0
                             ax.press = None
                             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, polar=self.polar)
                             ax.tick_params(labelsize=8)
                             ax.firsttime = True
                             ax.index = 0
                             ax.press = None
                             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%'):
                     '''
                     Add new axes to the given figure
                     '''
                     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):
                     '''
                     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 = numpy.median(deltas)
                     index = numpy.where(deltas > 5 * x_median)
                     if len(index[0]) != 0:
                         z_buffer[::, index[0], ::] = self.__missing
                         z_buffer = numpy.ma.masked_inside(z_buffer,
 .99 * self.__missing,
 .01 * self.__missing)
                     return x_buffer, y_buffer, z_buffer
                 def decimate(self):
                     # dx = int(len(self.x)/self.__MAXNUMX) + 1
                     dy = int(len(self.y) / self.decimation) + 1
                     # x = self.x[::dx]
                     x = self.x
                     y = self.y[::dy]
                     z = self.z[::, ::, ::dy]
                     return x, y, z
                 def format(self):
                     '''
                     Set min and max values, labels, ticks and titles
                     '''
                     for n, ax in enumerate(self.axes):
                         if ax.firsttime:
                             if self.xaxis != 'time':
                                 xmin = self.xmin
                                 xmax = self.xmax
                             else:
                                 xmin = self.tmin
                                 xmax = self.tmin + self.xrange*60*60
                                 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
                                 ax.xaxis.set_major_locator(LinearLocator(9))
                             ymin = self.ymin if self.ymin is not None else numpy.nanmin(self.y[numpy.isfinite(self.y)])
                             ymax = self.ymax if self.ymax is not None else numpy.nanmax(self.y[numpy.isfinite(self.y)])
                             ax.set_facecolor(self.bgcolor)
                             if self.xscale:
                                 ax.xaxis.set_major_formatter(FuncFormatter(
                                     lambda x, pos: '{0:g}'.format(x*self.xscale)))
                             if self.yscale:
                                 ax.yaxis.set_major_formatter(FuncFormatter(
                                     lambda x, pos: '{0:g}'.format(x*self.yscale)))
                             if self.xlabel is not None:
                                 ax.set_xlabel(self.xlabel)
                             if self.ylabel is not None:
                                 ax.set_ylabel(self.ylabel)
                             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:
                                 ax.cbar = plt.colorbar(
                                     ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
                                 ax.cbar.ax.tick_params(labelsize=8)
                                 ax.cbar.ax.press = None
                                 if self.cb_label:
                                     ax.cbar.set_label(self.cb_label, size=8)
                                 elif self.cb_labels:
                                     ax.cbar.set_label(self.cb_labels[n], size=8)
                             else:
                                 ax.cbar = None
                             ax.set_xlim(xmin, xmax)
                             ax.set_ylim(ymin, ymax)
                             ax.firsttime = False
                             if self.grid:
                                 ax.grid(True)
                         if not self.polar:
                             ax.set_title('{} {} {}'.format(
                                 self.titles[n],
                                 self.getDateTime(self.data.max_time).strftime(
                                     '%Y-%m-%d %H:%M:%S'),
                                 self.time_label),
                                 size=8)
                         else:
                             ax.set_title('{}'.format(self.titles[n]), size=8)
                             ax.set_ylim(0, 90)
                             ax.set_yticks(numpy.arange(0, 90, 20))
                             ax.yaxis.labelpad = 40
                     if self.firsttime:
                         for n, fig in enumerate(self.figures):
                             fig.subplots_adjust(**self.plots_adjust)
                         self.firsttime = False
                 def clear_figures(self):
                     '''
                     Reset axes for redraw plots
                     '''
                     for ax in self.axes+self.pf_axes+self.cb_axes:
                         ax.clear()
                         ax.firsttime = True
                         if hasattr(ax, 'cbar') and ax.cbar:
                             ax.cbar.remove()
                 def __plot(self):
                     '''
                     Main function to plot, format and save figures
                     '''
                     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)
                             fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
                             fig.canvas.manager.set_window_title(self.CODE)
                             continue
                         fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
                                                                              self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
                         fig.canvas.draw()
                         if self.show:
                             fig.show()
                             figpause(0.01)
                         if self.save:
                             self.save_figure(n)
                     if self.server:
                         self.send_to_server()
                 def __update(self, dataOut, timestamp):
                     '''
                     '''
                     metadata = {
                         'yrange': dataOut.heightList,
                         'interval': dataOut.timeInterval,
                         'channels': dataOut.channelList
                     }
                     data, meta = self.update(dataOut)
                     metadata.update(meta)
                     self.data.update(data, timestamp, metadata)
                 def save_figure(self, n):
                     '''
                     '''
                     if (self.data.max_time - self.save_time) <= self.save_period:
                         return
                     self.save_time = self.data.max_time
                     fig = self.figures[n]
                     if self.throttle == 0:
                         figname = os.path.join(
                             self.save,
                             self.save_code,
                             '{}_{}.png'.format(
                                 self.save_code,
                                 self.getDateTime(self.data.max_time).strftime(
                                     '%Y%m%d_%H%M%S'
                                     ),
                                 )
                             )
                         log.log('Saving figure: {}'.format(figname), self.name)
                         if not os.path.isdir(os.path.dirname(figname)):
                             os.makedirs(os.path.dirname(figname))
                         fig.savefig(figname)
                     figname = os.path.join(
                         self.save,
                         '{}_{}.png'.format(
                             self.save_code,
                             self.getDateTime(self.data.min_time).strftime(
                                 '%Y%m%d'
                                 ),
                             )
                         )
                     log.log('Saving figure: {}'.format(figname), self.name)
                     if not os.path.isdir(os.path.dirname(figname)):
                         os.makedirs(os.path.dirname(figname))
                     fig.savefig(figname)
                 def send_to_server(self):
                     '''
                     '''
                     if self.exp_code == None:
                         log.warning('Missing `exp_code` skipping sending to server...')
                     last_time = self.data.max_time
                     interval = last_time - self.sender_time
                     if interval < self.sender_period:
                         return
                     self.sender_time = last_time
                     attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
                     for attr in attrs:
                         value = getattr(self, attr)
                         if value:
                             if isinstance(value, (numpy.float32, numpy.float64)):
                                 value = round(float(value), 2)
                             self.data.meta[attr] = value
                     if self.colormap == 'jet':
                         self.data.meta['colormap'] = 'Jet'
                     elif 'RdBu' in self.colormap:
                         self.data.meta['colormap'] = 'RdBu'
                     else:
                         self.data.meta['colormap'] = 'Viridis'
                     self.data.meta['interval'] = int(interval)
                     self.sender_queue.append(last_time)
                     while True:
                         try:
                             tm = self.sender_queue.popleft()
                         except IndexError:
                             break
                         msg = self.data.jsonify(tm, self.save_code, self.plot_type)
                         self.socket.send_string(msg)
                         socks = dict(self.poll.poll(2000))
                         if socks.get(self.socket) == zmq.POLLIN:
                             reply = self.socket.recv_string()
                             if reply == 'ok':
                                 log.log("Response from server ok", self.name)
                                 time.sleep(0.1)
                                 continue
                             else:
                                 log.warning(
                                     "Malformed reply from server: {}".format(reply), self.name)
                         else:
                             log.warning(
                                 "No response from server, retrying...", self.name)
                         self.sender_queue.appendleft(tm)
                         self.socket.setsockopt(zmq.LINGER, 0)
                         self.socket.close()
                         self.poll.unregister(self.socket)
                         self.socket = self.context.socket(zmq.REQ)
                         self.socket.connect(self.server)
                         self.poll.register(self.socket, zmq.POLLIN)
                         break
                 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
                 def plot(self):
                     '''
                     Must be defined in the child class, the actual plotting method
                     '''
                     raise NotImplementedError
                 def update(self, dataOut):
                     '''
                     Must be defined in the child class, update self.data with new data
                     '''
                     data = {
                         self.CODE: getattr(dataOut, 'data_{}'.format(self.CODE))
                     }
                     meta = {}
                     return data, meta
                 def run(self, dataOut, **kwargs):
                     '''
                     Main plotting routine
                     '''
                     if self.isConfig is False:
                         self.__setup(**kwargs)
                         if self.localtime:
                             self.getDateTime = datetime.datetime.fromtimestamp
                         else:
                             self.getDateTime = datetime.datetime.utcfromtimestamp
                         self.data.setup()
                         self.isConfig = True
                         if self.server:
                             self.context = zmq.Context()
                             self.socket = self.context.socket(zmq.REQ)
                             self.socket.connect(self.server)
                             self.poll = zmq.Poller()
                             self.poll.register(self.socket, zmq.POLLIN)
                     tm = getattr(dataOut, self.attr_time)
                     if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60:
                         self.save_time = tm
                         self.__plot()
                         self.tmin += self.xrange*60*60
                         self.data.setup()
                         self.clear_figures()
                     self.__update(dataOut, tm)
                     if self.isPlotConfig is False:
                         self.__setup_plot()
                         self.isPlotConfig = True
                         if self.xaxis == 'time':
                             dt = self.getDateTime(tm)
                             if self.xmin is None:
                                 self.tmin = tm
                                 self.xmin = dt.hour
                             minutes = (self.xmin-int(self.xmin)) * 60
                             seconds = (minutes - int(minutes)) * 60
                             self.tmin = (dt.replace(hour=int(self.xmin), minute=int(minutes), second=int(seconds)) -
                                     datetime.datetime(1970, 1, 1)).total_seconds()
                             if self.localtime:
                                 self.tmin += time.timezone
                             if self.xmin is not None and self.xmax is not None:
                                 self.xrange = self.xmax - self.xmin
                     if self.throttle == 0:
                         self.__plot()
                     else:
                         self.__throttle_plot(self.__plot)#, coerce=coerce)
                 def close(self):
                     if self.data and not self.data.flagNoData:
                         self.save_time = 0
                         self.__plot()
                     if self.data and not self.data.flagNoData and self.pause:
                         figpause(10)

schainpy/model/graphics/jroplot_parameters.py +1 -2

             import os
             import datetime
             import numpy
             from schainpy.model.graphics.jroplot_base import Plot, plt
             from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot
             from schainpy.utils import log
             EARTH_RADIUS = 6.3710e3
             def ll2xy(lat1, lon1, lat2, lon2):
                 p = 0.017453292519943295
                 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
                     numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
                 r = 12742 * numpy.arcsin(numpy.sqrt(a))
                 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
                                       * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
                 theta = -theta + numpy.pi/2
                 return r*numpy.cos(theta), r*numpy.sin(theta)
             def km2deg(km):
                 '''
                 Convert distance in km to degrees
                 '''
                 return numpy.rad2deg(km/EARTH_RADIUS)
             class SpectralMomentsPlot(SpectraPlot):
                 '''
                 Plot for Spectral Moments
                 '''
                 CODE = 'spc_moments'
                 colormap = 'jet'
                 plot_type = 'pcolor'
             class SnrPlot(RTIPlot):
                 '''
                 Plot for SNR Data
                 '''
                 CODE = 'snr'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'snr': 10*numpy.log10(dataOut.data_snr)
                     }
                     return data, {}
             class DopplerPlot(RTIPlot):
                 '''
                 Plot for DOPPLER Data (1st moment)
                 '''
                 CODE = 'dop'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'dop': 10*numpy.log10(dataOut.data_dop)
                     }
                     return data, {}
             class PowerPlot(RTIPlot):
                 '''
                 Plot for Power Data (0 moment)
                 '''
                 CODE = 'pow'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'pow': 10*numpy.log10(dataOut.data_pow)
                     }
+                    print("data",data)
                     return data, {}
             class SpectralWidthPlot(RTIPlot):
                 '''
                 Plot for Spectral Width Data (2nd moment)
                 '''
                 CODE = 'width'
                 colormap = 'jet'
                 def update(self, dataOut):
                     data = {
                         'width': dataOut.data_width
                     }
                     return data, {}
             class SkyMapPlot(Plot):
                 '''
                 Plot for meteors detection data
                 '''
                 CODE = 'param'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.width = 7.2
                     self.height = 7.2
                     self.nplots = 1
                     self.xlabel = 'Zonal Zenith Angle (deg)'
                     self.ylabel = 'Meridional Zenith Angle (deg)'
                     self.polar = True
                     self.ymin = -180
                     self.ymax = 180
                     self.colorbar = False
                 def plot(self):
                     arrayParameters = numpy.concatenate(self.data['param'])
                     error = arrayParameters[:, -1]
                     indValid = numpy.where(error == 0)[0]
                     finalMeteor = arrayParameters[indValid, :]
                     finalAzimuth = finalMeteor[:, 3]
                     finalZenith = finalMeteor[:, 4]
                     x = finalAzimuth * numpy.pi / 180
                     y = finalZenith
                     ax = self.axes[0]
                     if ax.firsttime:
                         ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
                     else:
                         ax.plot.set_data(x, y)
                     dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
                     dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
                     title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
                                                                                                  dt2,
                                                                                                  len(x))
                     self.titles[0] = title
             class GenericRTIPlot(Plot):
                 '''
                 Plot for data_xxxx object
                 '''
                 CODE = 'param'
                 colormap = 'viridis'
                 plot_type = 'pcolorbuffer'
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     self.nrows = self.data.shape('param')[0]
                     self.nplots = self.nrows
                     self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95, 'top': 0.95})
                     if not self.xlabel:
                         self.xlabel = 'Time'
                     self.ylabel = 'Height [km]'
                     if not self.titles:
                         self.titles = ['Param {}'.format(x) for x in range(self.nrows)]
                 def update(self, dataOut):
                     data = {
                         'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
                     }
                     meta = {}
                     return data, meta
                 def plot(self):
                     # self.data.normalize_heights()
                     self.x = self.data.times
                     self.y = self.data.yrange
                     self.z = self.data['param']
                     self.z = numpy.ma.masked_invalid(self.z)
                     if self.decimation is None:
                         x, y, z = self.fill_gaps(self.x, self.y, self.z)
                     else:
                         x, y, z = self.fill_gaps(*self.decimate())
                     for n, ax in enumerate(self.axes):
                         self.zmax = self.zmax if self.zmax is not None else numpy.max(
                             self.z[n])
                         self.zmin = self.zmin if self.zmin is not None else numpy.min(
                             self.z[n])
                         if ax.firsttime:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             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]
                                                    )
             class PolarMapPlot(Plot):
                 '''
                 Plot for weather radar
                 '''
                 CODE = 'param'
                 colormap = 'seismic'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.width = 9
                     self.height = 8
                     self.mode = self.data.meta['mode']
                     if self.channels is not None:
                         self.nplots = len(self.channels)
                         self.nrows = len(self.channels)
                     else:
                         self.nplots = self.data.shape(self.CODE)[0]
                         self.nrows = self.nplots
                         self.channels = list(range(self.nplots))
                     if self.mode == 'E':
                         self.xlabel = 'Longitude'
                         self.ylabel = 'Latitude'
                     else:
                         self.xlabel = 'Range (km)'
                         self.ylabel = 'Height (km)'
                     self.bgcolor = 'white'
                     self.cb_labels = self.data.meta['units']
                     self.lat = self.data.meta['latitude']
                     self.lon = self.data.meta['longitude']
                     self.xmin, self.xmax = float(
                         km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
                     self.ymin, self.ymax = float(
                         km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
                     #  self.polar = True
                 def plot(self):
                     for n, ax in enumerate(self.axes):
                         data = self.data['param'][self.channels[n]]
                         zeniths = numpy.linspace(
 , self.data.meta['max_range'], data.shape[1])
                         if self.mode == 'E':
                             azimuths = -numpy.radians(self.data.yrange)+numpy.pi/2
                             r, theta = numpy.meshgrid(zeniths, azimuths)
                             x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
                                 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
                             x = km2deg(x) + self.lon
                             y = km2deg(y) + self.lat
                         else:
                             azimuths = numpy.radians(self.data.yrange)
                             r, theta = numpy.meshgrid(zeniths, azimuths)
                             x, y = r*numpy.cos(theta), r*numpy.sin(theta)
                         self.y = zeniths
                         if ax.firsttime:
                             if self.zlimits is not None:
                                 self.zmin, self.zmax = self.zlimits[n]
                             ax.plt = ax.pcolormesh(  # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 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(  # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
                                 vmin=self.zmin,
                                 vmax=self.zmax,
                                 cmap=self.cmaps[n])
                         if self.mode == 'A':
                             continue
                         # plot district names
                         f = open('/data/workspace/schain_scripts/distrito.csv')
                         for line in f:
                             label, lon, lat = [s.strip() for s in line.split(',') if s]
                             lat = float(lat)
                             lon = float(lon)
                             # ax.plot(lon, lat, '.b', ms=2)
                             ax.text(lon, lat, label.decode('utf8'), ha='center',
                                     va='bottom', size='8', color='black')
                         # plot limites
                         limites = []
                         tmp = []
                         for line in open('/data/workspace/schain_scripts/lima.csv'):
                             if '#' in line:
                                 if tmp:
                                     limites.append(tmp)
                                 tmp = []
                                 continue
                             values = line.strip().split(',')
                             tmp.append((float(values[0]), float(values[1])))
                         for points in limites:
                             ax.add_patch(
                                 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
                         # plot Cuencas
                         for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
                             f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
                             values = [line.strip().split(',') for line in f]
                             points = [(float(s[0]), float(s[1])) for s in values]
                             ax.add_patch(Polygon(points, ec='b', fc='none'))
                         # plot grid
                         for r in (15, 30, 45, 60):
                             ax.add_artist(plt.Circle((self.lon, self.lat),
                                                      km2deg(r), color='0.6', fill=False, lw=0.2))
                             ax.text(
                                 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
                                 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
                                 '{}km'.format(r),
                                 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
                     if self.mode == 'E':
                         title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
                         label = 'E{:02d}'.format(int(self.data.meta['elevation']))
                     else:
                         title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
                         label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
                     self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
                     self.titles = ['{} {}'.format(
                         self.data.parameters[x], title) for x in self.channels]

schainpy/model/graphics/jroplot_spectra.py +5 -4

             # Copyright (c) 2012-2020 Jicamarca Radio Observatory
             # All rights reserved.
             #
             # Distributed under the terms of the BSD 3-clause license.
             """Classes to plot Spectra data
             """
             import os
             import numpy
             from schainpy.model.graphics.jroplot_base import Plot, plt, log
             class SpectraPlot(Plot):
                 '''
                 Plot for Spectra data
                 '''
                 CODE = 'spc'
                 colormap = 'jet'
                 plot_type = 'pcolor'
                 buffering = False
-                channelList = None
+                channelList = []
                 def setup(self):
                     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.height = 2.6 * self.nrows
                     self.cb_label = 'dB'
                     if self.showprofile:
                         self.width = 4 * self.ncols
                     else:
                         self.width = 3.5 * self.ncols
                     self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
                     self.ylabel = 'Range [km]'
                 def update(self, dataOut):
                     if self.channelList == None:
                         self.channelList = dataOut.channelList
                     data = {}
                     meta = {}
                     spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
                     data['spc'] = spc
                     data['rti'] = dataOut.getPower()
                     data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
                     meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
                     if self.CODE == 'spc_moments':
                         data['moments'] = dataOut.moments
                     return data, meta
                 def plot(self):
                     if self.xaxis == "frequency":
                         x = self.data.xrange[0]
                         self.xlabel = "Frequency (kHz)"
                     elif self.xaxis == "time":
                         x = self.data.xrange[1]
                         self.xlabel = "Time (ms)"
                     else:
                         x = self.data.xrange[2]
                         self.xlabel = "Velocity (m/s)"
                     if self.CODE == 'spc_moments':
                         x = self.data.xrange[2]
                         self.xlabel = "Velocity (m/s)"
                     self.titles = []
                     y = self.data.yrange
                     self.y = y
                     data = self.data[-1]
                     z = data['spc']
                     for n, ax in enumerate(self.axes):
                         noise = data['noise'][n]
                         if self.CODE == 'spc_moments':
                             mean = data['moments'][n, 1]
                         if ax.firsttime:
                             self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
                             self.xmin = self.xmin if self.xmin else -self.xmax
                             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.plt_profile = self.pf_axes[n].plot(
                                     data['rti'][n], 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_moments':
                                 ax.plt_mean = ax.plot(mean, y, color='k')[0]
                         else:
                             ax.plt.set_array(z[n].T.ravel())
                             if self.showprofile:
                                 ax.plt_profile.set_data(data['rti'][n], y)
                                 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
                             if self.CODE == 'spc_moments':
                                 ax.plt_mean.set_data(mean, y)
                         self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
             class CrossSpectraPlot(Plot):
                 CODE = 'cspc'
                 colormap = 'jet'
                 plot_type = 'pcolor'
                 zmin_coh = None
                 zmax_coh = None
                 zmin_phase = None
                 zmax_phase = None
                 def setup(self):
                     self.ncols = 4
                     self.nplots = len(self.data.pairs) * 2
                     self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
                     self.width = 3.1 * self.ncols
                     self.height = 2.6 * self.nrows
                     self.ylabel = 'Range [km]'
                     self.showprofile = False
                     self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     spc = dataOut.data_spc
                     cspc = dataOut.data_cspc
                     meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
                     meta['pairs'] = dataOut.pairsList
                     tmp = []
                     for n, pair in enumerate(meta['pairs']):
                         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
                         tmp.append(coh)
                         tmp.append(phase)
                     data['cspc'] = numpy.array(tmp)
                     return data, meta
                 def plot(self):
                     if self.xaxis == "frequency":
                         x = self.data.xrange[0]
                         self.xlabel = "Frequency (kHz)"
                     elif self.xaxis == "time":
                         x = self.data.xrange[1]
                         self.xlabel = "Time (ms)"
                     else:
                         x = self.data.xrange[2]
                         self.xlabel = "Velocity (m/s)"
                     self.titles = []
                     y = self.data.yrange
                     self.y = y
                     data = self.data[-1]
                     cspc = data['cspc']
                     for n in range(len(self.data.pairs)):
                         pair = self.data.pairs[n]
                         coh = cspc[n*2]
                         phase = cspc[n*2+1]
                         ax = self.axes[2 * n]
                         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 = self.axes[2 * n + 1]
                         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]))
             class RTIPlot(Plot):
                 '''
                 Plot for RTI data
                 '''
                 CODE = 'rti'
                 colormap = 'jet'
                 plot_type = 'pcolorbuffer'
                 titles = None
-                channelList = None
+                channelList = []
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
+                    print("dataChannels ",self.data.channels)
                     self.nrows = len(self.data.channels)
                     self.nplots = len(self.data.channels)
                     self.ylabel = 'Range [km]'
                     self.xlabel = 'Time'
                     self.cb_label = 'dB'
                     self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
                     self.titles = ['{} Channel {}'.format(
                         self.CODE.upper(), x) for x in range(self.nplots)]
+                    print("SETUP")
                 def update(self, dataOut):
-                    if self.channelList == None:
+                    if len(self.channelList) == 0:
                         self.channelList = dataOut.channelList
                     data = {}
                     meta = {}
                     data['rti'] = dataOut.getPower()
                     data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
                     return data, meta
                 def plot(self):
                     self.x = self.data.times
                     self.y = self.data.yrange
                     self.z = self.data[self.CODE]
                     self.z = numpy.ma.masked_invalid(self.z)
                     if self.channelList != None:
                         self.titles = ['{} Channel {}'.format(
                             self.CODE.upper(), x) for x in self.channelList]
                     if self.decimation is None:
                         x, y, z = self.fill_gaps(self.x, self.y, self.z)
                     else:
                         x, y, z = self.fill_gaps(*self.decimate())
                     for n, ax in enumerate(self.axes):
                         self.zmin = self.zmin if self.zmin else numpy.min(self.z)
                         self.zmax = self.zmax if self.zmax else numpy.max(self.z)
                         data = self.data[-1]
                         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(
                                     data['rti'][n], self.y)[0]
                                 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], 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(data['rti'][n], self.y)
                                 ax.plot_noise.set_data(numpy.repeat(
                                     data['noise'][n], len(self.y)), self.y)
             class CoherencePlot(RTIPlot):
                 '''
                 Plot for Coherence data
                 '''
                 CODE = 'coh'
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     self.nrows = len(self.data.pairs)
                     self.nplots = len(self.data.pairs)
                     self.ylabel = 'Range [km]'
                     self.xlabel = 'Time'
                     self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
                     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.cb_label = 'Degrees'
                         self.titles = [
                             'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     data['coh'] = dataOut.getCoherence()
                     meta['pairs'] = dataOut.pairsList
                     return data, meta
             class PhasePlot(CoherencePlot):
                 '''
                 Plot for Phase map data
                 '''
                 CODE = 'phase'
                 colormap = 'seismic'
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     data['phase'] = dataOut.getCoherence(phase=True)
                     meta['pairs'] = dataOut.pairsList
                     return data, meta
             class NoisePlot(Plot):
                 '''
                 Plot for noise
                 '''
                 CODE = 'noise'
                 plot_type = 'scatterbuffer'
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     self.nrows = 1
                     self.nplots = 1
                     self.ylabel = 'Intensity [dB]'
                     self.xlabel = 'Time'
                     self.titles = ['Noise']
                     self.colorbar = False
                     self.plots_adjust.update({'right': 0.85 })
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
                     meta['yrange'] = numpy.array([])
                     return data, meta
                 def plot(self):
                     x = self.data.times
                     xmin = self.data.min_time
                     xmax = xmin + self.xrange * 60 * 60
                     Y = self.data['noise']
                     if self.axes[0].firsttime:
                         self.ymin = numpy.nanmin(Y) - 5
                         self.ymax = numpy.nanmax(Y) + 5
                         for ch in self.data.channels:
                             y = Y[ch]
                             self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
                         plt.legend(bbox_to_anchor=(1.18, 1.0))
                     else:
                         for ch in self.data.channels:
                             y = Y[ch]
                             self.axes[0].lines[ch].set_data(x, y)
             class PowerProfilePlot(Plot):
                 CODE = 'pow_profile'
                 plot_type = 'scatter'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.nplots = 1
                     self.height = 4
                     self.width = 3
                     self.ylabel = 'Range [km]'
                     self.xlabel = 'Intensity [dB]'
                     self.titles = ['Power Profile']
                     self.colorbar = False
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     data[self.CODE] = dataOut.getPower()
                     return data, meta
                 def plot(self):
                     y = self.data.yrange
                     self.y = y
                     x = self.data[-1][self.CODE]
                     if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
                     if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
                     if self.axes[0].firsttime:
                         for ch in self.data.channels:
                             self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
                         plt.legend()
                     else:
                         for ch in self.data.channels:
                             self.axes[0].lines[ch].set_data(x[ch], y)
             class SpectraCutPlot(Plot):
                 CODE = 'spc_cut'
                 plot_type = 'scatter'
                 buffering = False
                 def setup(self):
                     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 + 1.5
                     self.height = 3 * self.nrows
                     self.ylabel = 'Power [dB]'
                     self.colorbar = False
                     self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
                     data['spc'] = spc
                     meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
                     return data, meta
                 def plot(self):
                     if self.xaxis == "frequency":
                         x = self.data.xrange[0][1:]
                         self.xlabel = "Frequency (kHz)"
                     elif self.xaxis == "time":
                         x = self.data.xrange[1]
                         self.xlabel = "Time (ms)"
                     else:
                         x = self.data.xrange[2]
                         self.xlabel = "Velocity (m/s)"
                     self.titles = []
                     y = self.data.yrange
                     z = self.data[-1]['spc']
                     if self.height_index:
                         index = numpy.array(self.height_index)
                     else:
                         index = numpy.arange(0, len(y), int((len(y))/9))
                     for n, ax in enumerate(self.axes):
                         if ax.firsttime:
                             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 numpy.nanmin(z)
                             self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
                             ax.plt = ax.plot(x, z[n, :, index].T)
                             labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
                             self.figures[0].legend(ax.plt, labels, loc='center right')
                         else:
                             for i, line in enumerate(ax.plt):
                                 line.set_data(x, z[n, :, index[i]])
                         self.titles.append('CH {}'.format(n))
             class BeaconPhase(Plot):
                 __isConfig = None
                 __nsubplots = None
                 PREFIX = 'beacon_phase'
                 def __init__(self):
                     Plot.__init__(self)
                     self.timerange = 24*60*60
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.counter_imagwr = 0
                     self.WIDTH = 800
                     self.HEIGHT = 400
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.xdata = None
                     self.ydata = None
                     self.PLOT_CODE = BEACON_CODE
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                     self.filename_phase = None
                     self.figfile = None
                     self.xmin = None
                     self.xmax = None
                 def getSubplots(self):
                     ncol = 1
                     nrow = 1
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 7
                     colspan = 6
                     self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH+self.WIDTHPROF,
                                       heightplot = self.HEIGHT+self.HEIGHTPROF,
                                       show=show)
                     nrow, ncol = self.getSubplots()
                     self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
                 def save_phase(self, filename_phase):
                     f = open(filename_phase,'w+')
                     f.write('\n\n')
                     f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
                     f.write('DD MM YYYY  HH MM SS   pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
                     f.close()
                 def save_data(self, filename_phase, data, data_datetime):
                     f=open(filename_phase,'a')
                     timetuple_data = data_datetime.timetuple()
                     day = str(timetuple_data.tm_mday)
                     month = str(timetuple_data.tm_mon)
                     year = str(timetuple_data.tm_year)
                     hour = str(timetuple_data.tm_hour)
                     minute = str(timetuple_data.tm_min)
                     second = str(timetuple_data.tm_sec)
                     f.write(day+' '+month+' '+year+'  '+hour+' '+minute+' '+second+'   '+str(data[0])+'   '+str(data[1])+'   '+str(data[2])+'   '+str(data[3])+'\n')
                     f.close()
                 def plot(self):
                     log.warning('TODO: Not yet implemented...')
                 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
                         xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
                         timerange=None,
                         save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
                         server=None, folder=None, username=None, password=None,
                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
                     if dataOut.flagNoData:
                         return dataOut
                     if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
                         return
                     if pairsList == None:
                         pairsIndexList = dataOut.pairsIndexList[:10]
                     else:
                         pairsIndexList = []
                         for pair in pairsList:
                             if pair not in dataOut.pairsList:
                                 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
                             pairsIndexList.append(dataOut.pairsList.index(pair))
                     if pairsIndexList == []:
                         return
              #         if len(pairsIndexList) > 4:
              #             pairsIndexList = pairsIndexList[0:4]
                     hmin_index = None
                     hmax_index = None
                     if hmin != None and hmax != None:
                         indexes = numpy.arange(dataOut.nHeights)
                         hmin_list = indexes[dataOut.heightList >= hmin]
                         hmax_list = indexes[dataOut.heightList <= hmax]
                         if hmin_list.any():
                             hmin_index = hmin_list[0]
                         if hmax_list.any():
                             hmax_index = hmax_list[-1]+1
                     x = dataOut.getTimeRange()
                     thisDatetime = dataOut.datatime
                     title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     xlabel = "Local Time"
                     ylabel = "Phase (degrees)"
                     update_figfile = False
                     nplots = len(pairsIndexList)
                     #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
                     phase_beacon = numpy.zeros(len(pairsIndexList))
                     for i in range(nplots):
                         pair = dataOut.pairsList[pairsIndexList[i]]
                         ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
                         powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
                         powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
                         avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
                         phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
                         if dataOut.beacon_heiIndexList:
                             phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
                         else:
                             phase_beacon[i] = numpy.average(phase)
                     if not self.isConfig:
                         nplots = len(pairsIndexList)
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         if timerange != None:
                             self.timerange = timerange
                         self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
                         if ymin == None: ymin = 0
                         if ymax == None: ymax = 360
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.isConfig = True
                         self.figfile = figfile
                         self.xdata = numpy.array([])
                         self.ydata = numpy.array([])
                         update_figfile = True
                         #open file beacon phase
                         path = '%s%03d' %(self.PREFIX, self.id)
                         beacon_file = os.path.join(path,'%s.txt'%self.name)
                         self.filename_phase = os.path.join(figpath,beacon_file)
                         #self.save_phase(self.filename_phase)
                     #store data beacon phase
                     #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
                     self.setWinTitle(title)
                     title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
                     legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
                     axes = self.axesList[0]
                     self.xdata = numpy.hstack((self.xdata, x[0:1]))
                     if len(self.ydata)==0:
                         self.ydata = phase_beacon.reshape(-1,1)
                     else:
                         self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
                     axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
                                 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
                                 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
                                 XAxisAsTime=True, grid='both'
                                 )
                     self.draw()
                     if dataOut.ltctime >= self.xmax:
                         self.counter_imagwr = wr_period
                         self.isConfig = False
                         update_figfile = True
                     self.save(figpath=figpath,
                               figfile=figfile,
                               save=save,
                               ftp=ftp,
                               wr_period=wr_period,
                               thisDatetime=thisDatetime,
                               update_figfile=update_figfile)
                     return dataOut

schainpy/model/io/jroIO_base.py 0 0

NO CONTENT: modified file

schainpy/model/io/jroIO_kamisr.py +1 -3

             '''
             Created on Set 9, 2015
             @author: roj-idl71 Karim Kuyeng
             @update: 2021, Joab Apaza
             '''
             import os
             import sys
             import glob
             import fnmatch
             import datetime
             import time
             import re
             import h5py
             import numpy
             try:
                 from gevent import sleep
             except:
                 from time import sleep
             from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
             from schainpy.model.data.jrodata import Voltage
             from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
             from numpy import imag
             class AMISRReader(ProcessingUnit):
                 '''
                 classdocs
                 '''
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     ProcessingUnit.__init__(self)
                     self.set = None
                     self.subset = None
                     self.extension_file = '.h5'
                     self.dtc_str = 'dtc'
                     self.dtc_id = 0
                     self.status = True
                     self.isConfig = False
                     self.dirnameList = []
                     self.filenameList = []
                     self.fileIndex = None
                     self.flagNoMoreFiles = False
                     self.flagIsNewFile = 0
                     self.filename = ''
                     self.amisrFilePointer = None
                     self.realBeamCode = []
                     self.beamCodeMap = None
                     self.azimuthList = []
                     self.elevationList = []
                     self.dataShape = None
                     self.profileIndex = 0
                     self.beamCodeByFrame = None
                     self.radacTimeByFrame = None
                     self.dataset = None
                     self.__firstFile = True
                     self.buffer = None
                     self.timezone = 'ut'
                     self.__waitForNewFile = 20
                     self.__filename_online = None
                     #Is really necessary create the output object in the initializer
                     self.dataOut = Voltage()
                     self.dataOut.error=False
                 def setup(self,path=None,
                                 startDate=None,
                                 endDate=None,
                                 startTime=None,
                                 endTime=None,
                                 walk=True,
                                 timezone='ut',
                                 all=0,
                                 code = None,
                                 nCode = 0,
                                 nBaud = 0,
                                 online=False):
                     self.timezone = timezone
                     self.all = all
                     self.online = online
                     self.code = code
                     self.nCode = int(nCode)
                     self.nBaud = int(nBaud)
                     #self.findFiles()
                     if not(online):
                         #Busqueda de archivos offline
                         self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk)
                     else:
                         self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
                     if not(self.filenameList):
                         print("There is no files into the folder: %s"%(path))
                         sys.exit()
                     self.fileIndex = 0
                     self.readNextFile(online)
                     '''
                     Add code
                     '''
                     self.isConfig = True
                     # print("Setup Done")
                     pass
                 def readAMISRHeader(self,fp):
                     if self.isConfig and (not self.flagNoMoreFiles):
                         newShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
                         if self.dataShape != newShape and newShape != None:
                             print("\nNEW FILE HAS A DIFFERENT SHAPE")
                             print(self.dataShape,newShape,"\n")
                             return 0
                     else:
                         self.dataShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
                     header = 'Raw11/Data/RadacHeader'
                     self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
                     if (self.startDate> datetime.date(2021, 7, 15)): #Se cambió la forma de extracción de Apuntes el 17
                         self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode()
                         self.trueBeams = self.beamcodeFile.split("\n")
                         self.trueBeams.pop()#remove last
                         [self.realBeamCode.append(x) for x in self.trueBeams if x not in self.realBeamCode]
                         self.beamCode = [int(x, 16) for x in self.realBeamCode]
                     else:
                         _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes
                         self.beamCode = _beamCode[0,:]
                     if self.beamCodeMap == None:
                         self.beamCodeMap = fp['Setup/BeamcodeMap']
                         for beam in self.beamCode:
                             beamAziElev = numpy.where(self.beamCodeMap[:,0]==beam)
                             beamAziElev = beamAziElev[0].squeeze()
                             self.azimuthList.append(self.beamCodeMap[beamAziElev,1])
                             self.elevationList.append(self.beamCodeMap[beamAziElev,2])
                             #print("Beamssss: ",self.beamCodeMap[beamAziElev,1],self.beamCodeMap[beamAziElev,2])
                     #print(self.beamCode)
                     #self.code = fp.get(header+'/Code') # NOT USE FOR THIS
                     self.frameCount = fp.get(header+'/FrameCount')# NOT USE FOR THIS
                     self.modeGroup = fp.get(header+'/ModeGroup')# NOT USE FOR THIS
                     self.nsamplesPulse = fp.get(header+'/NSamplesPulse')# TO GET NSA OR USING DATA FOR THAT
                     self.pulseCount = fp.get(header+'/PulseCount')# NOT USE FOR THIS
                     self.radacTime = fp.get(header+'/RadacTime')# 1st TIME ON FILE ANDE CALCULATE THE REST WITH IPP*nindexprofile
                     self.timeCount = fp.get(header+'/TimeCount')# NOT USE FOR THIS
                     self.timeStatus = fp.get(header+'/TimeStatus')# NOT USE FOR THIS
                     self.rangeFromFile = fp.get('Raw11/Data/Samples/Range')
                     self.frequency =  fp.get('Rx/Frequency')
                     txAus = fp.get('Raw11/Data/Pulsewidth')
                     self.nblocks = self.pulseCount.shape[0] #nblocks
                     self.nprofiles = self.pulseCount.shape[1] #nprofile
                     self.nsa = self.nsamplesPulse[0,0] #ngates
                     self.nchannels = len(self.beamCode)
                     self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
                     #self.__waitForNewFile = self.nblocks  # wait depending on the number of blocks since each block is 1 sec
                     self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
                     #filling radar controller header parameters
                     self.__ippKm = self.ippSeconds *.15*1e6 # in km
                     self.__txA = (txAus.value)*.15 #(ipp[us]*.15km/1us) in km
                     self.__txB = 0
                     nWindows=1
                     self.__nSamples = self.nsa
                     self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
                     self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
                     #for now until understand why the code saved is different (code included even though code not in tuf file)
                     #self.__codeType = 0
                    # self.__nCode = None
                    # self.__nBaud = None
                     self.__code = self.code
                     self.__codeType = 0
                     if self.code != None:
                         self.__codeType = 1
                     self.__nCode = self.nCode
                     self.__nBaud = self.nBaud
                     #self.__code = 0
                     #filling system header parameters
                     self.__nSamples = self.nsa
                     self.newProfiles = self.nprofiles/self.nchannels
                     self.__channelList = list(range(self.nchannels))
                     self.__frequency = self.frequency[0][0]
                     return 1
                 def createBuffers(self):
                     pass
                 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
                     self.path = path
                     self.startDate = startDate
                     self.endDate = endDate
                     self.startTime = startTime
                     self.endTime = endTime
                     self.walk = walk
                 def __checkPath(self):
                     if os.path.exists(self.path):
                         self.status = 1
                     else:
                         self.status = 0
                         print('Path:%s does not exists'%self.path)
                     return
                 def __selDates(self, amisr_dirname_format):
                     try:
                         year = int(amisr_dirname_format[0:4])
                         month = int(amisr_dirname_format[4:6])
                         dom = int(amisr_dirname_format[6:8])
                         thisDate = datetime.date(year,month,dom)
                         if (thisDate>=self.startDate and thisDate <= self.endDate):
                             return amisr_dirname_format
                     except:
                         return None
                 def __findDataForDates(self,online=False):
                     if not(self.status):
                         return None
                     pat = '\d+.\d+'
                     dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
                     dirnameList = [x for x in dirnameList if x!=None]
                     dirnameList = [x.string for x in dirnameList]
                     if not(online):
                         dirnameList = [self.__selDates(x) for x in dirnameList]
                         dirnameList = [x for x in dirnameList if x!=None]
                     if len(dirnameList)>0:
                         self.status = 1
                         self.dirnameList = dirnameList
                         self.dirnameList.sort()
                     else:
                         self.status = 0
                         return None
                 def __getTimeFromData(self):
                     startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
                     endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
                     print('Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader))
                     print('........................................')
                     filter_filenameList = []
                     self.filenameList.sort()
                     #for i in range(len(self.filenameList)-1):
                     for i in range(len(self.filenameList)):
                         filename = self.filenameList[i]
                         fp = h5py.File(filename,'r')
                         time_str = fp.get('Time/RadacTimeString')
                         startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
                         #startDateTimeStr_File = "2019-12-16 09:21:11"
                         junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
                         startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
                         #endDateTimeStr_File = "2019-12-16 11:10:11"
                         endDateTimeStr_File = time_str[-1][-1].decode('UTF-8').split('.')[0]
                         junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
                         endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
                         fp.close()
                         #print("check time", startDateTime_File)
                         if self.timezone == 'lt':
                             startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
                             endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
                         if (endDateTime_File>=startDateTime_Reader and endDateTime_File<=endDateTime_Reader):
                             filter_filenameList.append(filename)
                         if (endDateTime_File>endDateTime_Reader):
                             break
                     filter_filenameList.sort()
                     self.filenameList = filter_filenameList
                     return 1
                 def __filterByGlob1(self, dirName):
                     filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
                     filter_files.sort()
                     filterDict = {}
                     filterDict.setdefault(dirName)
                     filterDict[dirName] = filter_files
                     return filterDict
                 def __getFilenameList(self, fileListInKeys, dirList):
                     for value in fileListInKeys:
                         dirName = list(value.keys())[0]
                         for file in value[dirName]:
                             filename = os.path.join(dirName, file)
                             self.filenameList.append(filename)
                 def __selectDataForTimes(self, online=False):
                     #aun no esta implementado el filtro for tiempo
                     if not(self.status):
                         return None
                     dirList = [os.path.join(self.path,x) for x in self.dirnameList]
                     fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
                     self.__getFilenameList(fileListInKeys, dirList)
                     if not(online):
                         #filtro por tiempo
                         if not(self.all):
                             self.__getTimeFromData()
                         if len(self.filenameList)>0:
                             self.status = 1
                             self.filenameList.sort()
                         else:
                             self.status = 0
                             return None
                     else:
                         #get the last file - 1
                         self.filenameList = [self.filenameList[-2]]
                     new_dirnameList = []
                     for dirname in self.dirnameList:
                         junk = numpy.array([dirname in x for x in self.filenameList])
                         junk_sum = junk.sum()
                         if junk_sum > 0:
                             new_dirnameList.append(dirname)
                     self.dirnameList = new_dirnameList
                     return 1
                 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
                                         endTime=datetime.time(23,59,59),walk=True):
                     if endDate ==None:
                      startDate = datetime.datetime.utcnow().date()
                      endDate = datetime.datetime.utcnow().date()
                     self.__setParameters(path=path, startDate=startDate, endDate=endDate,startTime = startTime,endTime=endTime, walk=walk)
                     self.__checkPath()
                     self.__findDataForDates(online=True)
                     self.dirnameList = [self.dirnameList[-1]]
                     self.__selectDataForTimes(online=True)
                     return
                 def searchFilesOffLine(self,
                                         path,
                                         startDate,
                                         endDate,
                                         startTime=datetime.time(0,0,0),
                                         endTime=datetime.time(23,59,59),
                                         walk=True):
                     self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
                     self.__checkPath()
                     self.__findDataForDates()
                     self.__selectDataForTimes()
                     for i in range(len(self.filenameList)):
                         print("%s" %(self.filenameList[i]))
                     return
                 def __setNextFileOffline(self):
                     try:
                         self.filename = self.filenameList[self.fileIndex]
                         self.amisrFilePointer = h5py.File(self.filename,'r')
                         self.fileIndex += 1
                     except:
                         self.flagNoMoreFiles = 1
                         print("No more Files")
                         return 0
                     self.flagIsNewFile = 1
                     print("Setting the file: %s"%self.filename)
                     return 1
                 def __setNextFileOnline(self):
                     filename = self.filenameList[0]
                     if self.__filename_online != None:
                         self.__selectDataForTimes(online=True)
                         filename = self.filenameList[0]
                         wait = 0
                         self.__waitForNewFile=300 ## DEBUG:
                         while self.__filename_online == filename:
                             print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
                             if wait == 5:
                                 self.flagNoMoreFiles = 1
                                 return 0
                             sleep(self.__waitForNewFile)
                             self.__selectDataForTimes(online=True)
                             filename = self.filenameList[0]
                             wait += 1
                     self.__filename_online = filename
                     self.amisrFilePointer = h5py.File(filename,'r')
                     self.flagIsNewFile = 1
                     self.filename = filename
                     print("Setting the file: %s"%self.filename)
                     return 1
                 def readData(self):
                     buffer = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
                     re = buffer[:,:,:,0]
                     im = buffer[:,:,:,1]
                     dataset = re + im*1j
                     self.radacTime = self.amisrFilePointer.get('Raw11/Data/RadacHeader/RadacTime')
                     timeset = self.radacTime[:,0]
                     return dataset,timeset
                 def reshapeData(self):
                 #self.beamCodeByPulse, self.beamCode, self.nblocks, self.nprofiles, self.nsa,
                     channels = self.beamCodeByPulse[0,:]
                     nchan = self.nchannels
                     #self.newProfiles = self.nprofiles/nchan #must be defined on filljroheader
                     nblocks = self.nblocks
                     nsamples = self.nsa
                     #Dimensions : nChannels, nProfiles, nSamples
                     new_block = numpy.empty((nblocks, nchan, numpy.int_(self.newProfiles), nsamples), dtype="complex64")
                     ############################################
                     for thisChannel in range(nchan):
                         new_block[:,thisChannel,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[thisChannel])[0],:]
                     new_block = numpy.transpose(new_block, (1,0,2,3))
                     new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
                     return new_block
                 def updateIndexes(self):
                     pass
                 def fillJROHeader(self):
                     #fill radar controller header
                     self.dataOut.radarControllerHeaderObj = RadarControllerHeader(ipp=self.__ippKm,
                                                                                   txA=self.__txA,
                                                                                   txB=0,
                                                                                   nWindows=1,
                                                                                   nHeights=self.__nSamples,
                                                                                   firstHeight=self.__firstHeight,
                                                                                   deltaHeight=self.__deltaHeight,
                                                                                   codeType=self.__codeType,
                                                                                   nCode=self.__nCode, nBaud=self.__nBaud,
                                                                                   code = self.__code,
                                                                                   fClock=1)
                     #fill system header
                     self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
                                                                 nProfiles=self.newProfiles,
                                                                 nChannels=len(self.__channelList),
                                                                 adcResolution=14,
                                                                 pciDioBusWidth=32)
                     self.dataOut.type = "Voltage"
                     self.dataOut.data = None
                     self.dataOut.dtype = numpy.dtype([('real','<i8'),('imag','<i8')])
             #        self.dataOut.nChannels = 0
             #        self.dataOut.nHeights = 0
                     self.dataOut.nProfiles = self.newProfiles*self.nblocks
                     #self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
                     ranges = numpy.reshape(self.rangeFromFile.value,(-1))
                     self.dataOut.heightList =  ranges/1000.0 #km
                     self.dataOut.channelList = self.__channelList
                     self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
             #        self.dataOut.channelIndexList = None
                     self.dataOut.azimuthList = numpy.array(self.azimuthList)
                     self.dataOut.elevationList = numpy.array(self.elevationList)
                     self.dataOut.codeList = numpy.array(self.beamCode)
                     #print(self.dataOut.elevationList)
                     self.dataOut.flagNoData = True
                     #Set to TRUE if the data is discontinuous
                     self.dataOut.flagDiscontinuousBlock = False
                     self.dataOut.utctime = None
                     #self.dataOut.timeZone = -5 #self.__timezone/60  #timezone like jroheader, difference in minutes between UTC and localtime
                     if self.timezone == 'lt':
                         self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
                     else:
                         self.dataOut.timeZone = 0 #by default time is UTC
                     self.dataOut.dstFlag = 0
                     self.dataOut.errorCount = 0
                     self.dataOut.nCohInt = 1
                     self.dataOut.flagDecodeData = False #asumo que la data esta decodificada
                     self.dataOut.flagDeflipData = False #asumo que la data esta sin flip
                     self.dataOut.flagShiftFFT = False
                     self.dataOut.ippSeconds = self.ippSeconds
                     #Time interval between profiles
                     #self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
                     self.dataOut.frequency = self.__frequency
                     self.dataOut.realtime = self.online
                     pass
                 def readNextFile(self,online=False):
                     if not(online):
                         newFile = self.__setNextFileOffline()
                     else:
                         newFile = self.__setNextFileOnline()
                     if not(newFile):
                         self.dataOut.error = True
                         return 0
                     if not self.readAMISRHeader(self.amisrFilePointer):
                         self.dataOut.error = True
                         return 0
                     self.createBuffers()
                     self.fillJROHeader()
                     #self.__firstFile = False
                     self.dataset,self.timeset = self.readData()
                     if self.endDate!=None:
                      endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
                      time_str = self.amisrFilePointer.get('Time/RadacTimeString')
                      startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
                      junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
                      startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
                      if self.timezone == 'lt':
                       startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
                      if (startDateTime_File>endDateTime_Reader):
                          return 0
                     self.jrodataset = self.reshapeData()
                     #----self.updateIndexes()
                     self.profileIndex = 0
                     return 1
                 def __hasNotDataInBuffer(self):
                     if self.profileIndex >= (self.newProfiles*self.nblocks):
                         return 1
                     return 0
                 def getData(self):
                     if self.flagNoMoreFiles:
                         self.dataOut.flagNoData = True
                         return 0
                     if self.__hasNotDataInBuffer():
                         if not (self.readNextFile(self.online)):
                             return 0
                     if self.dataset is None: # setear esta condicion cuando no hayan datos por leer
                         self.dataOut.flagNoData = True
                         return 0
                     #self.dataOut.data = numpy.reshape(self.jrodataset[self.profileIndex,:],(1,-1))
                     self.dataOut.data = self.jrodataset[:,self.profileIndex,:]
                     #print("R_t",self.timeset)
                     #self.dataOut.utctime = self.jrotimeset[self.profileIndex]
                     #verificar basic header de jro data y ver si es compatible con este valor
                     #self.dataOut.utctime = self.timeset + (self.profileIndex * self.ippSeconds * self.nchannels)
                     indexprof = numpy.mod(self.profileIndex, self.newProfiles)
                     indexblock = self.profileIndex/self.newProfiles
                     #print (indexblock, indexprof)
                     diffUTC = 1.8e4 #UTC diference from peru in seconds --Joab
                     diffUTC = 0
                     t_comp = (indexprof * self.ippSeconds * self.nchannels) + diffUTC #
                     #print("utc :",indexblock," __ ",t_comp)
                     #print(numpy.shape(self.timeset))
                     self.dataOut.utctime = self.timeset[numpy.int_(indexblock)] + t_comp
                     #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
                     #print(self.dataOut.utctime)
                     self.dataOut.profileIndex = self.profileIndex
                     #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
                     self.dataOut.flagNoData = False
             #         if indexprof == 0:
             #             print self.dataOut.utctime
                     self.profileIndex += 1
-                    #return self.dataOut.data
+                    return self.dataOut.data
                 def run(self, **kwargs):
                     '''
                     This method will be called many times so here you should put all your code
                     '''
                     #print("running kamisr")
                     if not self.isConfig:
                         self.setup(**kwargs)
                         self.isConfig = True
                     self.getData()
-                    #return(self.dataOut.data)
-                    return(self.dataOut)

schainpy/model/io/jroIO_param.py +27 -2

             import os
             import time
             import datetime
             import numpy
             import h5py
             import schainpy.admin
             from schainpy.model.data.jrodata import *
             from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
             from schainpy.model.io.jroIO_base import *
             from schainpy.utils import log
             class HDFReader(Reader, ProcessingUnit):
                 """Processing unit to read HDF5 format files
                 This unit reads HDF5 files created with `HDFWriter` operation contains
                 by default two groups Data and Metadata all variables would be saved as `dataOut`
                 attributes.
                 It is possible to read any HDF5 file by given the structure in the `description`
                 parameter, also you can add extra values to metadata with the parameter `extras`.
                 Parameters:
                 -----------
                 path : str
                     Path where files are located.
                 startDate : date
                     Start date of the files
                 endDate : list
                     End date of the files
                 startTime : time
                     Start time of the files
                 endTime : time
                     End time of the files
                 description : dict, optional
                     Dictionary with the description of the HDF5 file
                 extras : dict, optional
                     Dictionary with extra metadata to be be added to `dataOut`
                 Examples
                 --------
                 desc = {
                     'Data': {
                         'data_output': ['u', 'v', 'w'],
                         'utctime': 'timestamps',
                     }  ,
                     'Metadata': {
                         'heightList': 'heights'
                     }
                 }
                 desc = {
                     'Data': {
                         'data_output': 'winds',
                         'utctime': 'timestamps'
                     },
                     'Metadata': {
                         'heightList': 'heights'
                     }
                 }
                 extras = {
                     'timeZone': 300
                 }
                 reader = project.addReadUnit(
                     name='HDFReader',
                     path='/path/to/files',
                     startDate='2019/01/01',
                     endDate='2019/01/31',
                     startTime='00:00:00',
                     endTime='23:59:59',
                     # description=json.dumps(desc),
                     # extras=json.dumps(extras),
                     )
                 """
                 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
                 def __init__(self):
                     ProcessingUnit.__init__(self)
                     self.dataOut = Parameters()
                     self.ext = ".hdf5"
                     self.optchar = "D"
                     self.meta = {}
                     self.data = {}
                     self.open_file = h5py.File
                     self.open_mode = 'r'
                     self.description = {}
                     self.extras = {}
                     self.filefmt = "*%Y%j***"
                     self.folderfmt = "*%Y%j"
                     self.utcoffset = 0
                 def setup(self, **kwargs):
                     self.set_kwargs(**kwargs)
                     if not self.ext.startswith('.'):
                         self.ext = '.{}'.format(self.ext)
                     if self.online:
                         log.log("Searching files in online mode...", self.name)
                         for nTries in range(self.nTries):
                             fullpath = self.searchFilesOnLine(self.path, self.startDate,
                                 self.endDate, self.expLabel, self.ext, self.walk,
                                 self.filefmt, self.folderfmt)
+                            pathname, filename = os.path.split(fullpath)
+                            print(pathname,filename)
                             try:
                                 fullpath = next(fullpath)
                             except:
                                 fullpath = None
                             if fullpath:
                                 break
                             log.warning(
                                 'Waiting {} sec for a valid file in {}: try {} ...'.format(
                                     self.delay, self.path, nTries + 1),
                                 self.name)
                             time.sleep(self.delay)
                         if not(fullpath):
                             raise schainpy.admin.SchainError(
                                 'There isn\'t any valid file in {}'.format(self.path))
                         pathname, filename = os.path.split(fullpath)
                         self.year = int(filename[1:5])
                         self.doy = int(filename[5:8])
                         self.set = int(filename[8:11]) - 1
                     else:
                         log.log("Searching files in {}".format(self.path), self.name)
                         self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
                             self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
                     self.setNextFile()
                     return
                 def readFirstHeader(self):
                     '''Read metadata and data'''
                     self.__readMetadata()
                     self.__readData()
                     self.__setBlockList()
                     if 'type' in self.meta:
                         self.dataOut = eval(self.meta['type'])()
                     for attr in self.meta:
+                        print("attr: ", attr)
                         setattr(self.dataOut, attr, self.meta[attr])
                     self.blockIndex = 0
                     return
                 def __setBlockList(self):
                     '''
                     Selects the data within the times defined
                     self.fp
                     self.startTime
                     self.endTime
                     self.blockList
                     self.blocksPerFile
                     '''
                     startTime = self.startTime
                     endTime = self.endTime
                     thisUtcTime = self.data['utctime'] + self.utcoffset
                     self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
                     thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
+                    self.startFileDatetime = thisDatetime
+                    print("datee ",self.startFileDatetime)
                     thisDate = thisDatetime.date()
                     thisTime = thisDatetime.time()
                     startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
                     endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
                     ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
                     self.blockList = ind
                     self.blocksPerFile = len(ind)
+                    self.blocksPerFile = len(thisUtcTime)
                     return
                 def __readMetadata(self):
                     '''
                     Reads Metadata
                     '''
                     meta = {}
                     if self.description:
                         for key, value in self.description['Metadata'].items():
                             meta[key] = self.fp[value][()]
                     else:
                         grp = self.fp['Metadata']
                         for name in grp:
                             meta[name] = grp[name][()]
                     if self.extras:
                         for key, value in self.extras.items():
                             meta[key] = value
                     self.meta = meta
                     return
+                def checkForRealPath(self, nextFile, nextDay):
+                    # print("check FRP")
+                    # dt = self.startFileDatetime + datetime.timedelta(1)
+                    # filename = '{}.{}{}'.format(self.path, dt.strftime('%Y%m%d'), self.ext)
+                    # fullfilename = os.path.join(self.path, filename)
+                    # print("check Path ",fullfilename,filename)
+                    # if os.path.exists(fullfilename):
+                    #     return fullfilename, filename
+                    # return None, filename
+                    return None,None
                 def __readData(self):
                     data = {}
                     if self.description:
                         for key, value in self.description['Data'].items():
                             if isinstance(value, str):
                                 if isinstance(self.fp[value], h5py.Dataset):
                                     data[key] = self.fp[value][()]
                                 elif isinstance(self.fp[value], h5py.Group):
                                     array = []
                                     for ch in self.fp[value]:
                                         array.append(self.fp[value][ch][()])
                                     data[key] = numpy.array(array)
                             elif isinstance(value, list):
                                 array = []
                                 for ch in value:
                                     array.append(self.fp[ch][()])
                                 data[key] = numpy.array(array)
                     else:
                         grp = self.fp['Data']
                         for name in grp:
                             if isinstance(grp[name], h5py.Dataset):
                                 array = grp[name][()]
                             elif isinstance(grp[name], h5py.Group):
                                 array = []
                                 for ch in grp[name]:
                                     array.append(grp[name][ch][()])
                                 array = numpy.array(array)
                             else:
                                 log.warning('Unknown type: {}'.format(name))
                             if name in self.description:
                                 key = self.description[name]
                             else:
                                 key = name
                             data[key] = array
                     self.data = data
                     return
                 def getData(self):
+                    if not self.isDateTimeInRange(self.startFileDatetime, self.startDate, self.endDate, self.startTime, self.endTime):
+                        self.dataOut.flagNoData = True
+                        self.dataOut.error = True
+                        return
                     for attr in self.data:
                         if self.data[attr].ndim == 1:
                             setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
                         else:
                             setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
                     self.dataOut.flagNoData = False
                     self.blockIndex += 1
                     log.log("Block No. {}/{} -> {}".format(
                         self.blockIndex,
                         self.blocksPerFile,
                         self.dataOut.datatime.ctime()), self.name)
                     return
                 def run(self, **kwargs):
                     if not(self.isConfig):
                         self.setup(**kwargs)
                         self.isConfig = True
                     if self.blockIndex == self.blocksPerFile:
                         self.setNextFile()
                     self.getData()
                     return
             @MPDecorator
             class HDFWriter(Operation):
                 """Operation to write HDF5 files.
                 The HDF5 file contains by default two groups Data and Metadata where
                 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
                 parameters, data attributes are normaly time dependent where the metadata
                 are not.
                 It is possible to customize the structure of the HDF5 file with the
                 optional description parameter see the examples.
                 Parameters:
                 -----------
                 path : str
                     Path where files will be saved.
                 blocksPerFile : int
                     Number of blocks per file
                 metadataList : list
                     List of the dataOut attributes that will be saved as metadata
                 dataList : int
                     List of the dataOut attributes that will be saved as data
                 setType : bool
                     If True the name of the files corresponds to the timestamp of the data
                 description : dict, optional
                     Dictionary with the desired description of the HDF5 file
                 Examples
                 --------
                 desc = {
                     'data_output': {'winds': ['z', 'w', 'v']},
                     'utctime': 'timestamps',
                     'heightList': 'heights'
                 }
                 desc = {
                     'data_output': ['z', 'w', 'v'],
                     'utctime': 'timestamps',
                     'heightList': 'heights'
                 }
                 desc = {
                     'Data': {
                         'data_output': 'winds',
                         'utctime': 'timestamps'
                     },
                     'Metadata': {
                         'heightList': 'heights'
                     }
                 }
                 writer = proc_unit.addOperation(name='HDFWriter')
                 writer.addParameter(name='path', value='/path/to/file')
                 writer.addParameter(name='blocksPerFile', value='32')
                 writer.addParameter(name='metadataList', value='heightList,timeZone')
                 writer.addParameter(name='dataList',value='data_output,utctime')
                 # writer.addParameter(name='description',value=json.dumps(desc))
                 """
                 ext = ".hdf5"
                 optchar = "D"
                 filename = None
                 path = None
                 setFile = None
                 fp = None
                 firsttime = True
                 #Configurations
                 blocksPerFile = None
                 blockIndex = None
                 dataOut = None
                 #Data Arrays
                 dataList = None
                 metadataList = None
                 currentDay = None
                 lastTime = None
                 def __init__(self):
                     Operation.__init__(self)
                     return
                 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None):
                     self.path = path
                     self.blocksPerFile = blocksPerFile
                     self.metadataList = metadataList
                     self.dataList = [s.strip() for s in dataList]
                     self.setType = setType
                     self.description = description
                     if self.metadataList is None:
                         self.metadataList = self.dataOut.metadata_list
                     tableList = []
                     dsList = []
                     for i in range(len(self.dataList)):
                         dsDict = {}
                         if hasattr(self.dataOut, self.dataList[i]):
                             dataAux = getattr(self.dataOut, self.dataList[i])
                             dsDict['variable'] = self.dataList[i]
                         else:
                             log.warning('Attribute {} not found in dataOut', self.name)
                             continue
                         if dataAux is None:
                             continue
                         elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
                             dsDict['nDim'] = 0
                         else:
                             dsDict['nDim'] = len(dataAux.shape)
                             dsDict['shape'] = dataAux.shape
                             dsDict['dsNumber'] = dataAux.shape[0]
                             dsDict['dtype'] = dataAux.dtype
                         dsList.append(dsDict)
                     self.dsList = dsList
                     self.currentDay = self.dataOut.datatime.date()
                 def timeFlag(self):
                     currentTime = self.dataOut.utctime
                     timeTuple = time.localtime(currentTime)
                     dataDay = timeTuple.tm_yday
                     if self.lastTime is None:
                         self.lastTime = currentTime
                         self.currentDay = dataDay
                         return False
                     timeDiff = currentTime - self.lastTime
                     #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
                     if dataDay != self.currentDay:
                         self.currentDay = dataDay
                         return True
                     elif timeDiff > 3*60*60:
                         self.lastTime = currentTime
                         return True
                     else:
                         self.lastTime = currentTime
                         return False
                 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
                         dataList=[], setType=None, description={}):
                     self.dataOut = dataOut
                     if not(self.isConfig):
                         self.setup(path=path, blocksPerFile=blocksPerFile,
                                    metadataList=metadataList, dataList=dataList,
                                    setType=setType, description=description)
                         self.isConfig = True
                         self.setNextFile()
                     self.putData()
                     return
                 def setNextFile(self):
                     ext = self.ext
                     path = self.path
                     setFile = self.setFile
                     timeTuple = time.localtime(self.dataOut.utctime)
                     subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
                     fullpath = os.path.join(path, subfolder)
                     if os.path.exists(fullpath):
                         filesList = os.listdir(fullpath)
                         filesList = [k for k in filesList if k.startswith(self.optchar)]
                         if len( filesList ) > 0:
                             filesList = sorted(filesList, key=str.lower)
                             filen = filesList[-1]
                             # el filename debera tener el siguiente formato
                             # 0 1234 567 89A BCDE (hex)
                             # x YYYY DDD SSS .ext
                             if isNumber(filen[8:11]):
                                 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
                             else:
                                 setFile = -1
                         else:
                             setFile = -1 #inicializo mi contador de seteo
                     else:
                         os.makedirs(fullpath)
                         setFile = -1 #inicializo mi contador de seteo
                     if self.setType is None:
                         setFile += 1
                         file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
                                                        timeTuple.tm_year,
                                                        timeTuple.tm_yday,
                                                        setFile,
                                                        ext )
                     else:
                         setFile = timeTuple.tm_hour*60+timeTuple.tm_min
                         file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
                                                        timeTuple.tm_year,
                                                        timeTuple.tm_yday,
                                                        setFile,
                                                        ext )
                     self.filename = os.path.join( path, subfolder, file )
                     #Setting HDF5 File
                     self.fp = h5py.File(self.filename, 'w')
                     #write metadata
                     self.writeMetadata(self.fp)
                     #Write data
                     self.writeData(self.fp)
                 def getLabel(self, name, x=None):
                     if x is None:
                         if 'Data' in self.description:
                             data = self.description['Data']
                             if 'Metadata' in self.description:
                                 data.update(self.description['Metadata'])
                         else:
                             data = self.description
                         if name in data:
                             if isinstance(data[name], str):
                                 return data[name]
                             elif isinstance(data[name], list):
                                 return None
                             elif isinstance(data[name], dict):
                                 for key, value in data[name].items():
                                     return key
                         return name
                     else:
                         if 'Metadata' in self.description:
                             meta = self.description['Metadata']
                         else:
                             meta = self.description
                         if name in meta:
                             if isinstance(meta[name], list):
                                 return meta[name][x]
                             elif isinstance(meta[name], dict):
                                 for key, value in meta[name].items():
                                     return value[x]
                         if 'cspc' in name:
                             return 'pair{:02d}'.format(x)
                         else:
                             return 'channel{:02d}'.format(x)
                 def writeMetadata(self, fp):
                     if self.description:
                         if 'Metadata' in self.description:
                             grp = fp.create_group('Metadata')
                         else:
                             grp = fp
                     else:
                         grp = fp.create_group('Metadata')
                     for i in range(len(self.metadataList)):
                         if not hasattr(self.dataOut, self.metadataList[i]):
                             log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
                             continue
                         value = getattr(self.dataOut, self.metadataList[i])
                         if isinstance(value, bool):
                             if value is True:
                                 value = 1
                             else:
                                 value = 0
                         grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
                     return
                 def writeData(self, fp):
                     if self.description:
                         if 'Data' in self.description:
                             grp = fp.create_group('Data')
                         else:
                             grp = fp
                     else:
                         grp = fp.create_group('Data')
                     dtsets = []
                     data = []
                     for dsInfo in self.dsList:
                         if dsInfo['nDim'] == 0:
                             ds = grp.create_dataset(
                                 self.getLabel(dsInfo['variable']),
                                 (self.blocksPerFile, ),
                                 chunks=True,
                                 dtype=numpy.float64)
                             dtsets.append(ds)
                             data.append((dsInfo['variable'], -1))
                         else:
                             label = self.getLabel(dsInfo['variable'])
                             if label is not None:
                                 sgrp = grp.create_group(label)
                             else:
                                 sgrp = grp
                             for i in range(dsInfo['dsNumber']):
                                 ds = sgrp.create_dataset(
                                     self.getLabel(dsInfo['variable'], i),
                                     (self.blocksPerFile, ) + dsInfo['shape'][1:],
                                     chunks=True,
                                     dtype=dsInfo['dtype'])
                                 dtsets.append(ds)
                                 data.append((dsInfo['variable'], i))
                     fp.flush()
                     log.log('Creating file: {}'.format(fp.filename), self.name)
                     self.ds = dtsets
                     self.data = data
                     self.firsttime = True
                     self.blockIndex = 0
                     return
                 def putData(self):
                     if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
                         self.closeFile()
                         self.setNextFile()
                     for i, ds in enumerate(self.ds):
                         attr, ch = self.data[i]
                         if ch == -1:
                             ds[self.blockIndex] = getattr(self.dataOut, attr)
                         else:
                             ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
                     self.fp.flush()
                     self.blockIndex += 1
                     log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
                     return
                 def closeFile(self):
                     if self.blockIndex != self.blocksPerFile:
                         for ds in self.ds:
                             ds.resize(self.blockIndex, axis=0)
                     if self.fp:
                         self.fp.flush()
                         self.fp.close()
                 def close(self):
                     self.closeFile()

schainpy/model/proc/jroproc_parameters.py 0 0

	1	NO CONTENT: modified file			NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff

schainpy/model/proc/jroproc_spectra.py +31 -31

             # Copyright (c) 2012-2020 Jicamarca Radio Observatory
             # All rights reserved.
             #
             # Distributed under the terms of the BSD 3-clause license.
             """Spectra processing Unit and operations
             Here you will find the processing unit `SpectraProc` and several operations
             to work with Spectra data type
             """
             import time
             import itertools
             import numpy
             import math
             from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
             from schainpy.model.data.jrodata import Spectra
             from schainpy.model.data.jrodata import hildebrand_sekhon
             from schainpy.utils import log
             from scipy.optimize import curve_fit
             class SpectraProc(ProcessingUnit):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
                     self.buffer = None
                     self.firstdatatime = None
                     self.profIndex = 0
                     self.dataOut = Spectra()
                     self.id_min = None
                     self.id_max = None
                     self.setupReq = False #Agregar a todas las unidades de proc
                 def __updateSpecFromVoltage(self):
                     self.dataOut.timeZone = self.dataIn.timeZone
                     self.dataOut.dstFlag = self.dataIn.dstFlag
                     self.dataOut.errorCount = self.dataIn.errorCount
                     self.dataOut.useLocalTime = self.dataIn.useLocalTime
                     try:
                         self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
                     except:
                         pass
                     self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
                     self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
                     self.dataOut.channelList = self.dataIn.channelList
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
                     self.dataOut.nProfiles = self.dataOut.nFFTPoints
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
                     self.dataOut.utctime = self.firstdatatime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
                     self.dataOut.flagShiftFFT = False
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.nIncohInt = 1
                     self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
                     self.dataOut.frequency = self.dataIn.frequency
                     self.dataOut.realtime = self.dataIn.realtime
                     self.dataOut.azimuth = self.dataIn.azimuth
                     self.dataOut.zenith = self.dataIn.zenith
                     self.dataOut.codeList = self.dataIn.codeList
                     self.dataOut.azimuthList = self.dataIn.azimuthList
                     self.dataOut.elevationList = self.dataIn.elevationList
                 def __getFft(self):
                     """
                     Convierte valores de Voltaje a Spectra
                     Affected:
                         self.dataOut.data_spc
                         self.dataOut.data_cspc
                         self.dataOut.data_dc
                         self.dataOut.heightList
                         self.profIndex
                         self.buffer
                         self.dataOut.flagNoData
                     """
                     fft_volt = numpy.fft.fft(
                         self.buffer, n=self.dataOut.nFFTPoints, axis=1)
                     fft_volt = fft_volt.astype(numpy.dtype('complex'))
                     dc = fft_volt[:, 0, :]
                     # calculo de self-spectra
                     fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
                     spc = fft_volt * numpy.conjugate(fft_volt)
                     spc = spc.real
                     blocksize = 0
                     blocksize += dc.size
                     blocksize += spc.size
                     cspc = None
                     pairIndex = 0
                     if self.dataOut.pairsList != None:
                         # calculo de cross-spectra
                         cspc = numpy.zeros(
                             (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
                         for pair in self.dataOut.pairsList:
                             if pair[0] not in self.dataOut.channelList:
                                 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
                                     str(pair), str(self.dataOut.channelList)))
                             if pair[1] not in self.dataOut.channelList:
                                 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
                                     str(pair), str(self.dataOut.channelList)))
                             cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
                                 numpy.conjugate(fft_volt[pair[1], :, :])
                             pairIndex += 1
                         blocksize += cspc.size
                     self.dataOut.data_spc = spc
                     self.dataOut.data_cspc = cspc
                     self.dataOut.data_dc = dc
                     self.dataOut.blockSize = blocksize
                     self.dataOut.flagShiftFFT = False
                 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
                     if self.dataIn.type == "Spectra":
                         self.dataOut.copy(self.dataIn)
                         if shift_fft:
                             #desplaza a la derecha en el eje 2 determinadas posiciones
                             shift = int(self.dataOut.nFFTPoints/2)
                             self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
                             if self.dataOut.data_cspc is not None:
                                 #desplaza a la derecha en el eje 2 determinadas posiciones
                                 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
                         if pairsList:
                             self.__selectPairs(pairsList)
                     elif self.dataIn.type == "Voltage":
                         self.dataOut.flagNoData = True
                         if nFFTPoints == None:
                             raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
                         if nProfiles == None:
                             nProfiles = nFFTPoints
                         if ippFactor == None:
                             self.dataOut.ippFactor = 1
                         self.dataOut.nFFTPoints = nFFTPoints
                         if self.buffer is None:
                             self.buffer = numpy.zeros((self.dataIn.nChannels,
                                                        nProfiles,
                                                        self.dataIn.nHeights),
                                                       dtype='complex')
                         if self.dataIn.flagDataAsBlock:
                             nVoltProfiles = self.dataIn.data.shape[1]
                             if nVoltProfiles == nProfiles:
                                 self.buffer = self.dataIn.data.copy()
                                 self.profIndex = nVoltProfiles
                             elif nVoltProfiles < nProfiles:
                                 if self.profIndex == 0:
                                     self.id_min = 0
                                     self.id_max = nVoltProfiles
                                 self.buffer[:, self.id_min:self.id_max,
                                             :] = self.dataIn.data
                                 self.profIndex += nVoltProfiles
                                 self.id_min += nVoltProfiles
                                 self.id_max += nVoltProfiles
                             else:
                                 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
                                     self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
                                 self.dataOut.flagNoData = True
                         else:
                             self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
                             self.profIndex += 1
                         if self.firstdatatime == None:
                             self.firstdatatime = self.dataIn.utctime
                         if self.profIndex == nProfiles:
                             self.__updateSpecFromVoltage()
                             if pairsList == None:
                                 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
                             else:
                                 self.dataOut.pairsList = pairsList
                             self.__getFft()
                             self.dataOut.flagNoData = False
                             self.firstdatatime = None
                             self.profIndex = 0
                     else:
                         raise ValueError("The type of input object '%s' is not valid".format(
                             self.dataIn.type))
                 def __selectPairs(self, pairsList):
                     if not pairsList:
                         return
                     pairs = []
                     pairsIndex = []
                     for pair in pairsList:
                         if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
                             continue
                         pairs.append(pair)
                         pairsIndex.append(pairs.index(pair))
                     self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
                     self.dataOut.pairsList = pairs
                     return
                 def selectFFTs(self, minFFT, maxFFT ):
                     """
                     Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
                     minFFT<= FFT <= maxFFT
                     """
                     if (minFFT > maxFFT):
                         raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
                     if (minFFT < self.dataOut.getFreqRange()[0]):
                         minFFT = self.dataOut.getFreqRange()[0]
                     if (maxFFT > self.dataOut.getFreqRange()[-1]):
                         maxFFT = self.dataOut.getFreqRange()[-1]
                     minIndex = 0
                     maxIndex = 0
                     FFTs = self.dataOut.getFreqRange()
                     inda = numpy.where(FFTs >= minFFT)
                     indb = numpy.where(FFTs <= maxFFT)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(FFTs)
                     self.selectFFTsByIndex(minIndex, maxIndex)
                     return 1
                 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
                     newheis = numpy.where(
                         self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
                     if hei_ref != None:
                         newheis = numpy.where(self.dataOut.heightList > hei_ref)
                     minIndex = min(newheis[0])
                     maxIndex = max(newheis[0])
                     data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
                     heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
                     # determina indices
                     nheis = int(self.dataOut.radarControllerHeaderObj.txB /
                                 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
                     avg_dB = 10 * \
                         numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
                     beacon_dB = numpy.sort(avg_dB)[-nheis:]
                     beacon_heiIndexList = []
                     for val in avg_dB.tolist():
                         if val >= beacon_dB[0]:
                             beacon_heiIndexList.append(avg_dB.tolist().index(val))
                     #data_spc = data_spc[:,:,beacon_heiIndexList]
                     data_cspc = None
                     if self.dataOut.data_cspc is not None:
                         data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
                         #data_cspc = data_cspc[:,:,beacon_heiIndexList]
                     data_dc = None
                     if self.dataOut.data_dc is not None:
                         data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
                         #data_dc = data_dc[:,beacon_heiIndexList]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_cspc = data_cspc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.heightList = heightList
                     self.dataOut.beacon_heiIndexList = beacon_heiIndexList
                     return 1
                 def selectFFTsByIndex(self, minIndex, maxIndex):
                     """
                     """
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
                     if (maxIndex >= self.dataOut.nProfiles):
                         maxIndex = self.dataOut.nProfiles-1
                     #Spectra
                     data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
                     data_cspc = None
                     if self.dataOut.data_cspc is not None:
                         data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
                     data_dc = None
                     if self.dataOut.data_dc is not None:
                         data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_cspc = data_cspc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
                     self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
                     self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
                     return 1
                 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
                     # validacion de rango
                     if minHei == None:
                         minHei = self.dataOut.heightList[0]
                     if maxHei == None:
                         maxHei = self.dataOut.heightList[-1]
                     if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
                         print('minHei: %.2f is out of the heights range' % (minHei))
                         print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
                         minHei = self.dataOut.heightList[0]
                     if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
                         print('maxHei: %.2f is out of the heights range' % (maxHei))
                         print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
                         maxHei = self.dataOut.heightList[-1]
                     # validacion de velocidades
                     velrange = self.dataOut.getVelRange(1)
                     if minVel == None:
                         minVel = velrange[0]
                     if maxVel == None:
                         maxVel = velrange[-1]
                     if (minVel < velrange[0]) or (minVel > maxVel):
                         print('minVel: %.2f is out of the velocity range' % (minVel))
                         print('minVel is setting to %.2f' % (velrange[0]))
                         minVel = velrange[0]
                     if (maxVel > velrange[-1]) or (maxVel < minVel):
                         print('maxVel: %.2f is out of the velocity range' % (maxVel))
                         print('maxVel is setting to %.2f' % (velrange[-1]))
                         maxVel = velrange[-1]
                     # seleccion de indices para rango
                     minIndex = 0
                     maxIndex = 0
                     heights = self.dataOut.heightList
                     inda = numpy.where(heights >= minHei)
                     indb = numpy.where(heights <= maxHei)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(heights)
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError("some value in (%d,%d) is not valid" % (
                             minIndex, maxIndex))
                     if (maxIndex >= self.dataOut.nHeights):
                         maxIndex = self.dataOut.nHeights - 1
                     # seleccion de indices para velocidades
                     indminvel = numpy.where(velrange >= minVel)
                     indmaxvel = numpy.where(velrange <= maxVel)
                     try:
                         minIndexVel = indminvel[0][0]
                     except:
                         minIndexVel = 0
                     try:
                         maxIndexVel = indmaxvel[0][-1]
                     except:
                         maxIndexVel = len(velrange)
                     # seleccion del espectro
                     data_spc = self.dataOut.data_spc[:,
                                                      minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
                     # estimacion de ruido
                     noise = numpy.zeros(self.dataOut.nChannels)
                     for channel in range(self.dataOut.nChannels):
                         daux = data_spc[channel, :, :]
                         sortdata = numpy.sort(daux, axis=None)
                         noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
                     self.dataOut.noise_estimation = noise.copy()
                     return 1
             class removeDC(Operation):
                 def run(self, dataOut, mode=2):
                     self.dataOut = dataOut
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     if jcspectra is not None:
                         jcspectraExist = True
                         num_pairs = jcspectra.shape[0]
                     else:
                         jcspectraExist = False
                     freq_dc = int(jspectra.shape[1] / 2)
                     ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
                     ind_vel = ind_vel.astype(int)
                     if ind_vel[0] < 0:
                         ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
                     if mode == 1:
                         jspectra[:, freq_dc, :] = (
                             jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2  # CORRECCION
                         if jcspectraExist:
                             jcspectra[:, freq_dc, :] = (
                                 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
                     if mode == 2:
                         vel = numpy.array([-2, -1, 1, 2])
                         xx = numpy.zeros([4, 4])
                         for fil in range(4):
                             xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
                         xx_inv = numpy.linalg.inv(xx)
                         xx_aux = xx_inv[0, :]
                         for ich in range(num_chan):
                             yy = jspectra[ich, ind_vel, :]
                             jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
                             junkid = jspectra[ich, freq_dc, :] <= 0
                             cjunkid = sum(junkid)
                             if cjunkid.any():
                                 jspectra[ich, freq_dc, junkid.nonzero()] = (
                                     jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
                         if jcspectraExist:
                             for ip in range(num_pairs):
                                 yy = jcspectra[ip, ind_vel, :]
                                 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
                     self.dataOut.data_spc = jspectra
                     self.dataOut.data_cspc = jcspectra
                     return self.dataOut
             # import matplotlib.pyplot as plt
             def fit_func( x, a0, a1, a2): #, a3, a4, a5):
                 z = (x - a1) / a2
                 y = a0 * numpy.exp(-z**2 / a2)  #+ a3 + a4 * x + a5 * x**2
                 return y
             class CleanRayleigh(Operation):
                 def __init__(self):
                     Operation.__init__(self)
                     self.i=0
                     self.isConfig = False
                     self.__dataReady = False
                     self.__profIndex = 0
                     self.byTime = False
                     self.byProfiles = False
                     self.bloques =  None
                     self.bloque0 = None
                     self.index = 0
                     self.buffer = 0
                     self.buffer2 = 0
                     self.buffer3 = 0
                 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
                     self.nChannels = dataOut.nChannels
                     self.nProf = dataOut.nProfiles
                     self.nPairs = dataOut.data_cspc.shape[0]
                     self.pairsArray = numpy.array(dataOut.pairsList)
                     self.spectra = dataOut.data_spc
                     self.cspectra = dataOut.data_cspc
                     self.heights = dataOut.heightList #alturas totales
                     self.nHeights =  len(self.heights)
                     self.min_hei = min_hei
                     self.max_hei = max_hei
                     if (self.min_hei == None):
                         self.min_hei = 0
                     if (self.max_hei == None):
                         self.max_hei = dataOut.heightList[-1]
                     self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
                     self.heightsClean = self.heights[self.hval] #alturas filtradas
                     self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
                     self.nHeightsClean = len(self.heightsClean)
                     self.channels = dataOut.channelList
                     self.nChan = len(self.channels)
                     self.nIncohInt = dataOut.nIncohInt
                     self.__initime = dataOut.utctime
                     self.maxAltInd = self.hval[-1]+1
                     self.minAltInd = self.hval[0]
                     self.crosspairs = dataOut.pairsList
                     self.nPairs = len(self.crosspairs)
                     self.normFactor = dataOut.normFactor
                     self.nFFTPoints = dataOut.nFFTPoints
                     self.ippSeconds = dataOut.ippSeconds
                     self.currentTime = self.__initime
                     self.pairsArray = numpy.array(dataOut.pairsList)
                     self.factor_stdv = factor_stdv
-                    print("CHANNELS: ",[x for x in self.channels])
+                    #print("CHANNELS: ",[x for x in self.channels])
                     if n != None :
                         self.byProfiles = True
                         self.nIntProfiles = n
                     else:
                         self.__integrationtime = timeInterval
                     self.__dataReady = False
                     self.isConfig = True
                 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
                     #print (dataOut.utctime)
                     if not self.isConfig :
                         #print("Setting config")
                         self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
                         #print("Config Done")
                     tini=dataOut.utctime
                     if self.byProfiles:
                         if self.__profIndex ==  self.nIntProfiles:
                             self.__dataReady = True
                     else:
                         if (tini - self.__initime) >= self.__integrationtime:
                             #print(tini - self.__initime,self.__profIndex)
                             self.__dataReady = True
                             self.__initime = tini
                     #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
                     if self.__dataReady:
-                        print("Data ready",self.__profIndex)
+                        #print("Data ready",self.__profIndex)
                         self.__profIndex = 0
                         jspc = self.buffer
                         jcspc = self.buffer2
                         #jnoise = self.buffer3
                         self.buffer = dataOut.data_spc
                         self.buffer2 = dataOut.data_cspc
                         #self.buffer3 = dataOut.noise
                         self.currentTime = dataOut.utctime
                         if numpy.any(jspc) :
                             #print( jspc.shape, jcspc.shape)
                             jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
                             jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
                             self.__dataReady = False
                             #print( jspc.shape, jcspc.shape)
                             dataOut.flagNoData = False
                         else:
                             dataOut.flagNoData = True
                             self.__dataReady = False
                             return dataOut
                     else:
                         #print( len(self.buffer))
                         if numpy.any(self.buffer):
                             self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
                             self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
                             self.buffer3 += dataOut.data_dc
                         else:
                             self.buffer = dataOut.data_spc
                             self.buffer2 = dataOut.data_cspc
                             self.buffer3 = dataOut.data_dc
                         #print self.index, self.fint
                         #print self.buffer2.shape
                         dataOut.flagNoData = True                   ## NOTE: ?? revisar LUEGO
                         self.__profIndex += 1
                         return dataOut                             ## NOTE: REV
                     #index = tini.tm_hour*12+tini.tm_min/5
                     '''REVISAR'''
                     # jspc = jspc/self.nFFTPoints/self.normFactor
                     # jcspc = jcspc/self.nFFTPoints/self.normFactor
-                    #dataOut.data_spc,dataOut.data_cspc = self.CleanRayleigh(dataOut,jspc,jcspc,crosspairs,heights,channels,nProf,nHei,nChan,nPairs,nIncohInt,nBlocks=nBlocks)
-                    #tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.min_hei,self.max_hei)
-                    tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
-                    #jspectra = tmp_spectra*len(jspc[:,0,0,0])
-                    #jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
+                    tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
                     dataOut.data_spc = tmp_spectra
                     dataOut.data_cspc = tmp_cspectra
+                    #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
                     dataOut.data_dc = self.buffer3
                     dataOut.nIncohInt *= self.nIntProfiles
                     dataOut.utctime = self.currentTime             #tiempo promediado
                     #print("Time: ",time.localtime(dataOut.utctime))
                     # dataOut.data_spc = sat_spectra
                     # dataOut.data_cspc = sat_cspectra
                     self.buffer = 0
                     self.buffer2 = 0
                     self.buffer3 = 0
                     return dataOut
                 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
-                    print("OP cleanRayleigh")
+                    #print("OP cleanRayleigh")
                     #import matplotlib.pyplot as plt
                     #for k in range(149):
                     rfunc = cspectra.copy() #self.bloques
-                    val_spc = spectra*0.0 #self.bloque0*0.0
+                    #rfunc = cspectra
-                    val_cspc = cspectra*0.0 #self.bloques*0.0
+                    #val_spc = spectra*0.0 #self.bloque0*0.0
-                    in_sat_spectra = spectra.copy()  #self.bloque0
+                    #val_cspc = cspectra*0.0 #self.bloques*0.0
-                    in_sat_cspectra = cspectra.copy()  #self.bloques
+                    #in_sat_spectra = spectra.copy()  #self.bloque0
+                    #in_sat_cspectra = cspectra.copy()  #self.bloques
-                    raxs = math.ceil(math.sqrt(self.nPairs))
+                    #raxs = math.ceil(math.sqrt(self.nPairs))
-                    caxs = math.ceil(self.nPairs/raxs)
+                    #caxs = math.ceil(self.nPairs/raxs)
                     #print(self.hval)
                     #print numpy.absolute(rfunc[:,0,0,14])
+                    gauss_fit, covariance = None, None
                     for ih in range(self.minAltInd,self.maxAltInd):
                         for ifreq in range(self.nFFTPoints):
                             # fig, axs = plt.subplots(raxs, caxs)
                             # fig2, axs2 = plt.subplots(raxs, caxs)
-                            col_ax = 0
+                            # col_ax = 0
-                            row_ax = 0
+                            # row_ax = 0
+                            #print(len(self.nPairs))
                             for ii in range(self.nPairs):                       #PARES DE CANALES SELF y CROSS
                                 #print("ii: ",ii)
-                                if (col_ax%caxs==0 and col_ax!=0):
+                                # if (col_ax%caxs==0 and col_ax!=0):
-                                    col_ax = 0
+                                #     col_ax = 0
-                                    row_ax += 1
+                                #     row_ax += 1
                                 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))           #Potencia?
                                 #print(func2clean.shape)
                                 val = (numpy.isfinite(func2clean)==True).nonzero()
                                 if len(val)>0:                                                              #limitador
                                     min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
                                     if min_val <= -40 :
                                         min_val = -40
                                     max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
                                     if max_val >= 200 :
                                         max_val = 200
                                     #print min_val, max_val
                                     step = 1
                                     #print("Getting bins and the histogram")
                                     x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
                                     y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
                                     #print(len(y_dist),len(binstep[:-1]))
                                     #print(row_ax,col_ax, " ..")
                                     #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
                                     mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
                                     sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
                                     parg = [numpy.amax(y_dist),mean,sigma]
-                                    gauss_fit, covariance = None, None
-                                    newY = None
+                                    #newY = None
                                     try :
                                         gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
                                         mode = gauss_fit[1]
                                         stdv = gauss_fit[2]
                                         #print(" FIT OK",gauss_fit)
                                         '''
                                         newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
                                         axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
                                         axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
                                         axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))'''
                                     except:
                                         mode = mean
                                         stdv = sigma
                                         #print("FIT FAIL")
                                     #print(mode,stdv)
-                                    #Removing echoes greater than mode + 3*stdv
+                                    #Removing echoes greater than mode + std_factor*stdv
-                                    #factor_stdv = 2
                                     noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
                                     #noval tiene los indices que se van a remover
                                     #print("Pair ",ii," novals: ",len(noval[0]))
                                     if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
                                         novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
                                         #print(novall)
                                         #print("   ",self.pairsArray[ii])
                                         cross_pairs = self.pairsArray[ii]
                                             #Getting coherent echoes which are removed.
                                         # if len(novall[0]) > 0:
                                         #
                                         #     val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
                                         #     val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
                                         #     val_cspc[novall[0],ii,ifreq,ih] = 1
                                             #print("OUT NOVALL 1")
                                             #Removing coherent from ISR data
                                         chA = self.channels.index(cross_pairs[0])
                                         chB = self.channels.index(cross_pairs[1])
                                         new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0])
-                                        mean_cspc = numpy.mean(new_a)
+                                        cspectra[noval,ii,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
                                         new_b = numpy.delete(spectra[:,chA,ifreq,ih], noval[0])
-                                        mean_spc0 = numpy.mean(new_b)
+                                        spectra[noval,chA,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
                                         new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
-                                        mean_spc1 = numpy.mean(new_c)
+                                        spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
-                                        spectra[noval,chA,ifreq,ih] = mean_spc0
-                                        spectra[noval,chB,ifreq,ih] = mean_spc1
-                                        cspectra[noval,ii,ifreq,ih] = mean_cspc
                                     '''
                                     func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
                                     y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
                                     axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
                                     axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
                                     axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
                                     '''
-                                col_ax += 1    #contador de ploteo columnas
+                                #col_ax += 1    #contador de ploteo columnas
                                 ##print(col_ax)
                                 '''
                             title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+"   Alt: "+str(self.heights[ih])+ " km"
                             title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+"   Alt: "+str(self.heights[ih])+ " km  CLEANED"
                             fig.suptitle(title)
                             fig2.suptitle(title2)
                             plt.show()'''
                         '''                channels = channels
                                         cross_pairs = cross_pairs
                                             #print("OUT NOVALL 2")
                                         vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
                                         vcross1 = (cross_pairs[1] == channels[ii]).nonzero()
                                         vcross = numpy.concatenate((vcross0,vcross1),axis=None)
                                             #print('vcros =', vcross)
                                             #Getting coherent echoes which are removed.
                                         if len(novall) > 0:
                                                 #val_spc[novall,ii,ifreq,ih] = 1
                                              val_spc[ii,ifreq,ih,novall] = 1
                                              if len(vcross) > 0:
                                                 val_cspc[vcross,ifreq,ih,novall] = 1
                                             #Removing coherent from ISR data.
                                         self.bloque0[ii,ifreq,ih,noval] = numpy.nan
                                         if len(vcross) > 0:
                                             self.bloques[vcross,ifreq,ih,noval] = numpy.nan
                         '''
-                    print("Getting average of the spectra and cross-spectra from incoherent echoes.")
+                    #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
                     out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
                     out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
                     for ih in range(self.nHeights):
                         for ifreq in range(self.nFFTPoints):
                             for ich in range(self.nChan):
                                 tmp = spectra[:,ich,ifreq,ih]
                                 valid = (numpy.isfinite(tmp[:])==True).nonzero()
             #                     if ich == 0 and ifreq == 0 and ih == 17 :
             #                         print tmp
             #                         print valid
             #                         print len(valid[0])
                                 #print('TMP',tmp)
                                 if len(valid[0]) >0 :
                                     out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
                                 #for icr in range(nPairs):
                             for icr in range(self.nPairs):
                                 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
                                 valid = (numpy.isfinite(tmp)==True).nonzero()
                                 if len(valid[0]) > 0:
                                     out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
                     '''
                    # print('##########################################################')
                     print("Removing fake coherent echoes (at least 4 points around the point)")
                     val_spectra = numpy.sum(val_spc,0)
                     val_cspectra = numpy.sum(val_cspc,0)
                     val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
                     val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
                     for i in range(nChan):
                         for j in range(nProf):
                             for k in range(nHeights):
                                 if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 :
                                     val_spc[:,i,j,k] = 0.0
                     for i in range(nPairs):
                         for j in range(nProf):
                             for k in range(nHeights):
                                 if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
                                     val_cspc[:,i,j,k] = 0.0
                     # val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProf*nHeights*nChan))
                     # if numpy.isfinite(val_spectra)==str(True):
                     #    noval = (val_spectra<1).nonzero()
                     #    if len(noval) > 0:
                     #        val_spc[:,noval] = 0.0
                     #        val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHeights))
                         #val_cspc = numpy.reshape(val_spc, (149,nChan*nHeights*nProf))
                         #if numpy.isfinite(val_cspectra)==str(True):
                          #   noval = (val_cspectra<1).nonzero()
                           #  if len(noval) > 0:
                            #     val_cspc[:,noval] = 0.0
                             #    val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHeights))
                     tmp_sat_spectra = spectra.copy()
                     tmp_sat_spectra = tmp_sat_spectra*numpy.nan
                     tmp_sat_cspectra = cspectra.copy()
                     tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
                     '''
                     # fig = plt.figure(figsize=(6,5))
                     # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
                     # ax = fig.add_axes([left, bottom, width, height])
                     # cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:])))
                     # ax.clabel(cp, inline=True,fontsize=10)
                     # plt.show()
                     '''
                     val = (val_spc > 0).nonzero()
                     if len(val[0]) > 0:
                         tmp_sat_spectra[val] = in_sat_spectra[val]
                     val = (val_cspc > 0).nonzero()
                     if len(val[0]) > 0:
                         tmp_sat_cspectra[val] = in_sat_cspectra[val]
                     print("Getting average of the spectra and cross-spectra from incoherent echoes 2")
                     sat_spectra = numpy.zeros((nChan,nProf,nHeights), dtype=float)
                     sat_cspectra = numpy.zeros((nPairs,nProf,nHeights), dtype=complex)
                     for ih in range(nHeights):
                         for ifreq in range(nProf):
                             for ich in range(nChan):
                                 tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
                                 valid = (numpy.isfinite(tmp)).nonzero()
                                 if len(valid[0]) > 0:
                                     sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
                             for icr in range(nPairs):
                                 tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih])
                                 valid = (numpy.isfinite(tmp)).nonzero()
                                 if len(valid[0]) > 0:
                                     sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
                     '''
                         #self.__dataReady= True
                         #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra
                     #if not self.__dataReady:
                         #return None, None
                     #return out_spectra, out_cspectra ,sat_spectra,sat_cspectra
                     return out_spectra, out_cspectra
                 def REM_ISOLATED_POINTS(self,array,rth):
             #         import matplotlib.pyplot as plt
                     if rth == None :
                         rth = 4
                     print("REM ISO")
                     num_prof = len(array[0,:,0])
                     num_hei = len(array[0,0,:])
                     n2d = len(array[:,0,0])
                     for ii in range(n2d) :
                         #print ii,n2d
                         tmp = array[ii,:,:]
                         #print tmp.shape, array[ii,101,:],array[ii,102,:]
                       # fig = plt.figure(figsize=(6,5))
                       # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
                       # ax = fig.add_axes([left, bottom, width, height])
                       # x = range(num_prof)
                       # y = range(num_hei)
                       # cp = ax.contour(y,x,tmp)
                       # ax.clabel(cp, inline=True,fontsize=10)
                       # plt.show()
                         #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
                         tmp = numpy.reshape(tmp,num_prof*num_hei)
                         indxs1 = (numpy.isfinite(tmp)==True).nonzero()
                         indxs2 = (tmp > 0).nonzero()
                         indxs1 = (indxs1[0])
                         indxs2 = indxs2[0]
                         #indxs1 = numpy.array(indxs1[0])
                         #indxs2 = numpy.array(indxs2[0])
                         indxs = None
                         #print indxs1 , indxs2
                         for iv in range(len(indxs2)):
                             indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
                             #print len(indxs2), indv
                             if len(indv[0]) > 0  :
                                 indxs =  numpy.concatenate((indxs,indxs2[iv]), axis=None)
               #           print indxs
                         indxs = indxs[1:]
                         #print(indxs, len(indxs))
                         if len(indxs) < 4 :
                             array[ii,:,:] = 0.
                             return
                         xpos = numpy.mod(indxs ,num_hei)
                         ypos = (indxs / num_hei)
                         sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
                         #print sx
                         xpos = xpos[sx]
                         ypos = ypos[sx]
                  # *********************************** Cleaning isolated points **********************************
                         ic = 0
                         while True :
                             r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
                             #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
                             #plt.plot(r)
                             #plt.show()
                             no_coh1 = (numpy.isfinite(r)==True).nonzero()
                             no_coh2 = (r <= rth).nonzero()
                             #print r, no_coh1, no_coh2
                             no_coh1 = numpy.array(no_coh1[0])
                             no_coh2 = numpy.array(no_coh2[0])
                             no_coh = None
                             #print valid1 , valid2
                             for iv in range(len(no_coh2)):
                                 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
                                 if len(indv[0]) > 0  :
                                     no_coh =  numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
                             no_coh = no_coh[1:]
                             #print len(no_coh), no_coh
                             if len(no_coh) < 4 :
                                 #print xpos[ic], ypos[ic], ic
                 #                plt.plot(r)
                 #                plt.show()
                                 xpos[ic] = numpy.nan
                                 ypos[ic] = numpy.nan
                             ic = ic + 1
                             if  (ic == len(indxs)) :
                                 break
                             #print( xpos, ypos)
                         indxs = (numpy.isfinite(list(xpos))==True).nonzero()
                         #print indxs[0]
                         if len(indxs[0]) < 4 :
                             array[ii,:,:] = 0.
                             return
                         xpos = xpos[indxs[0]]
                         ypos = ypos[indxs[0]]
                         for i in range(0,len(ypos)):
                             ypos[i]=int(ypos[i])
                         junk = tmp
                         tmp = junk*0.0
                         tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
                         array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
                         #print array.shape
                         #tmp = numpy.reshape(tmp,(num_prof,num_hei))
                         #print tmp.shape
                         # fig = plt.figure(figsize=(6,5))
                         # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
                         # ax = fig.add_axes([left, bottom, width, height])
                         # x = range(num_prof)
                         # y = range(num_hei)
                         # cp = ax.contour(y,x,array[ii,:,:])
                         # ax.clabel(cp, inline=True,fontsize=10)
                         # plt.show()
                     return array
             class removeInterference(Operation):
                 def removeInterference2(self):
                     cspc = self.dataOut.data_cspc
                     spc = self.dataOut.data_spc
                     Heights = numpy.arange(cspc.shape[2])
                     realCspc = numpy.abs(cspc)
                     for i in range(cspc.shape[0]):
                         LinePower= numpy.sum(realCspc[i], axis=0)
                         Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
                         SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
                         InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
                         InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
                         InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
                         InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
                         #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
                         if len(InterferenceRange)<int(cspc.shape[1]*0.3):
                             cspc[i,InterferenceRange,:] = numpy.NaN
                     self.dataOut.data_cspc = cspc
                 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     jnoise = self.dataOut.getNoise()
                     num_incoh = self.dataOut.nIncohInt
                     num_channel = jspectra.shape[0]
                     num_prof = jspectra.shape[1]
                     num_hei = jspectra.shape[2]
                     # hei_interf
                     if hei_interf is None:
                         count_hei = int(num_hei / 2)
                         hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
                         hei_interf = numpy.asarray(hei_interf)[0]
                     # nhei_interf
                     if (nhei_interf == None):
                         nhei_interf = 5
                     if (nhei_interf < 1):
                         nhei_interf = 1
                     if (nhei_interf > count_hei):
                         nhei_interf = count_hei
                     if (offhei_interf == None):
                         offhei_interf = 0
                     ind_hei = list(range(num_hei))
             #         mask_prof = numpy.asarray(range(num_prof - 2)) + 1
             #         mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
                     mask_prof = numpy.asarray(list(range(num_prof)))
                     num_mask_prof = mask_prof.size
                     comp_mask_prof = [0, num_prof / 2]
                     # noise_exist:    Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
                     if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
                         jnoise = numpy.nan
                     noise_exist = jnoise[0] < numpy.Inf
                     # Subrutina de Remocion de la Interferencia
                     for ich in range(num_channel):
                         # Se ordena los espectros segun su potencia (menor a mayor)
                         power = jspectra[ich, mask_prof, :]
                         power = power[:, hei_interf]
                         power = power.sum(axis=0)
                         psort = power.ravel().argsort()
                         # Se estima la interferencia promedio en los Espectros de Potencia empleando
                         junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
                             offhei_interf, nhei_interf + offhei_interf))]]]
                         if noise_exist:
                             #    tmp_noise = jnoise[ich] / num_prof
                             tmp_noise = jnoise[ich]
                         junkspc_interf = junkspc_interf - tmp_noise
                         #junkspc_interf[:,comp_mask_prof] = 0
                         jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
                         jspc_interf = jspc_interf.transpose()
                         # Calculando el espectro de interferencia promedio
                         noiseid = numpy.where(
                             jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
                         noiseid = noiseid[0]
                         cnoiseid = noiseid.size
                         interfid = numpy.where(
                             jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
                         interfid = interfid[0]
                         cinterfid = interfid.size
                         if (cnoiseid > 0):
                             jspc_interf[noiseid] = 0
                         # Expandiendo los perfiles a limpiar
                         if (cinterfid > 0):
                             new_interfid = (
                                 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
                             new_interfid = numpy.asarray(new_interfid)
                             new_interfid = {x for x in new_interfid}
                             new_interfid = numpy.array(list(new_interfid))
                             new_cinterfid = new_interfid.size
                         else:
                             new_cinterfid = 0
                         for ip in range(new_cinterfid):
                             ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
                             jspc_interf[new_interfid[ip]
                                         ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
                         jspectra[ich, :, ind_hei] = jspectra[ich, :,
                                                              ind_hei] - jspc_interf  # Corregir indices
                         # Removiendo la interferencia del punto de mayor interferencia
                         ListAux = jspc_interf[mask_prof].tolist()
                         maxid = ListAux.index(max(ListAux))
                         if cinterfid > 0:
                             for ip in range(cinterfid * (interf == 2) - 1):
                                 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
                                        (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
                                 cind = len(ind)
                                 if (cind > 0):
                                     jspectra[ich, interfid[ip], ind] = tmp_noise * \
                                         (1 + (numpy.random.uniform(cind) - 0.5) /
                                          numpy.sqrt(num_incoh))
                             ind = numpy.array([-2, -1, 1, 2])
                             xx = numpy.zeros([4, 4])
                             for id1 in range(4):
                                 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
                             xx_inv = numpy.linalg.inv(xx)
                             xx = xx_inv[:, 0]
                             ind = (ind + maxid + num_mask_prof) % num_mask_prof
                             yy = jspectra[ich, mask_prof[ind], :]
                             jspectra[ich, mask_prof[maxid], :] = numpy.dot(
                                 yy.transpose(), xx)
                         indAux = (jspectra[ich, :, :] < tmp_noise *
                                   (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
                         jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
                             (1 - 1 / numpy.sqrt(num_incoh))
                     # Remocion de Interferencia en el Cross Spectra
                     if jcspectra is None:
                         return jspectra, jcspectra
                     num_pairs = int(jcspectra.size / (num_prof * num_hei))
                     jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
                     for ip in range(num_pairs):
                         #-------------------------------------------
                         cspower = numpy.abs(jcspectra[ip, mask_prof, :])
                         cspower = cspower[:, hei_interf]
                         cspower = cspower.sum(axis=0)
                         cspsort = cspower.ravel().argsort()
                         junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
                             offhei_interf, nhei_interf + offhei_interf))]]]
                         junkcspc_interf = junkcspc_interf.transpose()
                         jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
                         ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
                         median_real = int(numpy.median(numpy.real(
                             junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
                         median_imag = int(numpy.median(numpy.imag(
                             junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
                         comp_mask_prof = [int(e) for e in comp_mask_prof]
                         junkcspc_interf[comp_mask_prof, :] = numpy.complex(
                             median_real, median_imag)
                         for iprof in range(num_prof):
                             ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
                             jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
                         # Removiendo la Interferencia
                         jcspectra[ip, :, ind_hei] = jcspectra[ip,
                                                               :, ind_hei] - jcspc_interf
                         ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
                         maxid = ListAux.index(max(ListAux))
                         ind = numpy.array([-2, -1, 1, 2])
                         xx = numpy.zeros([4, 4])
                         for id1 in range(4):
                             xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
                         xx_inv = numpy.linalg.inv(xx)
                         xx = xx_inv[:, 0]
                         ind = (ind + maxid + num_mask_prof) % num_mask_prof
                         yy = jcspectra[ip, mask_prof[ind], :]
                         jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
                     # Guardar Resultados
                     self.dataOut.data_spc = jspectra
                     self.dataOut.data_cspc = jcspectra
                     return 1
                 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
                     self.dataOut = dataOut
                     if mode == 1:
                         self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
                     elif mode == 2:
                         self.removeInterference2()
                     return self.dataOut
             class IncohInt(Operation):
                 __profIndex = 0
                 __withOverapping = False
                 __byTime = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer_spc = None
                 __buffer_cspc = None
                 __buffer_dc = None
                 __dataReady = False
                 __timeInterval = None
                 n = None
                 def __init__(self):
                     Operation.__init__(self)
                 def setup(self, n=None, timeInterval=None, overlapping=False):
                     """
                     Set the parameters of the integration class.
                     Inputs:
                         n        :    Number of coherent integrations
                         timeInterval   :    Time of integration. If the parameter "n" is selected this one does not work
                         overlapping    :
                     """
                     self.__initime = None
                     self.__lastdatatime = 0
                     self.__buffer_spc = 0
                     self.__buffer_cspc = 0
                     self.__buffer_dc = 0
                     self.__profIndex = 0
                     self.__dataReady = False
                     self.__byTime = False
                     if n is None and timeInterval is None:
                         raise ValueError("n or timeInterval should be specified ...")
                     if n is not None:
                         self.n = int(n)
                     else:
                         self.__integrationtime = int(timeInterval)
                         self.n = None
                         self.__byTime = True
                 def putData(self, data_spc, data_cspc, data_dc):
                     """
                     Add a profile to the __buffer_spc and increase in one the __profileIndex
                     """
                     self.__buffer_spc += data_spc
                     if data_cspc is None:
                         self.__buffer_cspc = None
                     else:
                         self.__buffer_cspc += data_cspc
                     if data_dc is None:
                         self.__buffer_dc = None
                     else:
                         self.__buffer_dc += data_dc
                     self.__profIndex += 1
                     return
                 def pushData(self):
                     """
                     Return the sum of the last profiles and the profiles used in the sum.
                     Affected:
                     self.__profileIndex
                     """
                     data_spc = self.__buffer_spc
                     data_cspc = self.__buffer_cspc
                     data_dc = self.__buffer_dc
                     n = self.__profIndex
                     self.__buffer_spc = 0
                     self.__buffer_cspc = 0
                     self.__buffer_dc = 0
                     self.__profIndex = 0
                     return data_spc, data_cspc, data_dc, n
                 def byProfiles(self, *args):
                     self.__dataReady = False
                     avgdata_spc = None
                     avgdata_cspc = None
                     avgdata_dc = None
                     self.putData(*args)
                     if self.__profIndex == self.n:
                         avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
                         self.n = n
                         self.__dataReady = True
                     return avgdata_spc, avgdata_cspc, avgdata_dc
                 def byTime(self, datatime, *args):
                     self.__dataReady = False
                     avgdata_spc = None
                     avgdata_cspc = None
                     avgdata_dc = None
                     self.putData(*args)
                     if (datatime - self.__initime) >= self.__integrationtime:
                         avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
                         self.n = n
                         self.__dataReady = True
                     return avgdata_spc, avgdata_cspc, avgdata_dc
                 def integrate(self, datatime, *args):
                     if self.__profIndex == 0:
                         self.__initime = datatime
                     if self.__byTime:
                         avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
                             datatime, *args)
                     else:
                         avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
                     if not self.__dataReady:
                         return None, None, None, None
                     return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
                 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
                     if n == 1:
                         return dataOut
                     dataOut.flagNoData = True
                     if not self.isConfig:
                         self.setup(n, timeInterval, overlapping)
                         self.isConfig = True
                     avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
                                                                                         dataOut.data_spc,
                                                                                         dataOut.data_cspc,
                                                                                         dataOut.data_dc)
                     if self.__dataReady:
                         dataOut.data_spc = avgdata_spc
                         dataOut.data_cspc = avgdata_cspc
                         dataOut.data_dc = avgdata_dc
                         dataOut.nIncohInt *= self.n
                         dataOut.utctime = avgdatatime
                         dataOut.flagNoData = False
                     return dataOut
             class dopplerFlip(Operation):
                 def run(self, dataOut):
                     # arreglo 1: (num_chan, num_profiles, num_heights)
                     self.dataOut = dataOut
                     # JULIA-oblicua, indice 2
                     # arreglo 2: (num_profiles, num_heights)
                     jspectra = self.dataOut.data_spc[2]
                     jspectra_tmp = numpy.zeros(jspectra.shape)
                     num_profiles = jspectra.shape[0]
                     freq_dc = int(num_profiles / 2)
                     # Flip con for
                     for j in range(num_profiles):
                         jspectra_tmp[num_profiles-j-1]= jspectra[j]
                     # Intercambio perfil de DC con perfil inmediato anterior
                     jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
                     jspectra_tmp[freq_dc]= jspectra[freq_dc]
                     # canal modificado es re-escrito en el arreglo de canales
                     self.dataOut.data_spc[2] = jspectra_tmp
                     return self.dataOut

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages