schain Commit - r1653:f661d02a9309 · Jicamarca Repository

jroplot_parameters.py jroplot_spectra.py jroIO_digitalRF.py jroproc_parameters.py sophy_proc.py

eynilupu -

r1653:f661d02a9309

parent child

Context file:

r1653:f661d02a9309

Collapse all files

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

schainpy/model/graphics/jroplot_parameters.py +20 -9

             import os
             import datetime
             import warnings
             import numpy
             from mpl_toolkits.axisartist.grid_finder import FixedLocator, DictFormatter
             from matplotlib.patches import Circle
             import cartopy.crs as ccrs
             from cartopy.feature import ShapelyFeature
             import cartopy.io.shapereader as shpreader
             from schainpy.model.graphics.jroplot_base import Plot, plt
             from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot, SpectraCutPlot
             from schainpy.utils import log
             from schainpy.model.graphics.plotting_codes import cb_tables
             EARTH_RADIUS = 6.3710e3
             def antenna_to_cartesian(ranges, azimuths, elevations):
                 """
                 Return Cartesian coordinates from antenna coordinates.
                 Parameters
                 ----------
                 ranges : array
                     Distances to the center of the radar gates (bins) in kilometers.
                 azimuths : array
                     Azimuth angle of the radar in degrees.
                 elevations : array
                     Elevation angle of the radar in degrees.
                 Returns
                 -------
                 x, y, z : array
                     Cartesian coordinates in meters from the radar.
                 Notes
                 -----
                 The calculation for Cartesian coordinate is adapted from equations
 .28(b) and 2.28(c) of Doviak and Zrnic [1]_ assuming a
                 standard atmosphere (4/3 Earth's radius model).
                 .. math::
                     z = \\sqrt{r^2+R^2+2*r*R*sin(\\theta_e)} - R
                     s = R * arcsin(\\frac{r*cos(\\theta_e)}{R+z})
                     x = s * sin(\\theta_a)
                     y = s * cos(\\theta_a)
                 Where r is the distance from the radar to the center of the gate,
                 :math:`\\theta_a` is the azimuth angle, :math:`\\theta_e` is the
                 elevation angle, s is the arc length, and R is the effective radius
                 of the earth, taken to be 4/3 the mean radius of earth (6371 km).
                 References
                 ----------
                 .. [1] Doviak and Zrnic, Doppler Radar and Weather Observations, Second
                     Edition, 1993, p. 21.
                 """
                 theta_e = numpy.deg2rad(elevations) # elevation angle in radians.
                 theta_a = numpy.deg2rad(azimuths) # azimuth angle in radians.
                 R = 6371.0 * 1000.0 * 4.0 / 3.0     # effective radius of earth in meters.
                 r = ranges * 1000.0                 # distances to gates in meters.
                 z = (r ** 2 + R ** 2 + 2.0 * r * R * numpy.sin(theta_e)) ** 0.5 - R
                 s = R * numpy.arcsin(r * numpy.cos(theta_e) / (R + z))  # arc length in m.
                 x = s * numpy.sin(theta_a)
                 y = s * numpy.cos(theta_a)
                 return x, y, z
             def cartesian_to_geographic_aeqd(x, y, lon_0, lat_0, R=EARTH_RADIUS):
                 """
                 Azimuthal equidistant Cartesian to geographic coordinate transform.
                 Transform a set of Cartesian/Cartographic coordinates (x, y) to
                 geographic coordinate system (lat, lon) using a azimuthal equidistant
                 map projection [1]_.
                 .. math::
                     lat = \\arcsin(\\cos(c) * \\sin(lat_0) +
                                    (y * \\sin(c) * \\cos(lat_0) / \\rho))
                     lon = lon_0 + \\arctan2(
                         x * \\sin(c),
                         \\rho * \\cos(lat_0) * \\cos(c) - y * \\sin(lat_0) * \\sin(c))
                     \\rho = \\sqrt(x^2 + y^2)
                     c = \\rho / R
                 Where x, y are the Cartesian position from the center of projection;
                 lat, lon the corresponding latitude and longitude; lat_0, lon_0 are the
                 latitude and longitude of the center of the projection; R is the radius of
                 the earth (defaults to ~6371 km). lon is adjusted to be between -180 and
 .
                 Parameters
                 ----------
                 x, y : array-like
                     Cartesian coordinates in the same units as R, typically meters.
                 lon_0, lat_0 : float
                     Longitude and latitude, in degrees, of the center of the projection.
                 R : float, optional
                     Earth radius in the same units as x and y. The default value is in
                     units of meters.
                 Returns
                 -------
                 lon, lat : array
                     Longitude and latitude of Cartesian coordinates in degrees.
                 References
                 ----------
                 .. [1] Snyder, J. P. Map Projections--A Working Manual. U. S. Geological
                     Survey Professional Paper 1395, 1987, pp. 191-202.
                 """
                 x = numpy.atleast_1d(numpy.asarray(x))
                 y = numpy.atleast_1d(numpy.asarray(y))
                 lat_0_rad = numpy.deg2rad(lat_0)
                 lon_0_rad = numpy.deg2rad(lon_0)
                 rho = numpy.sqrt(x*x + y*y)
                 c = rho / R
                 with warnings.catch_warnings():
                     # division by zero may occur here but is properly addressed below so
                     # the warnings can be ignored
                     warnings.simplefilter("ignore", RuntimeWarning)
                     lat_rad = numpy.arcsin(numpy.cos(c) * numpy.sin(lat_0_rad) +
                                         y * numpy.sin(c) * numpy.cos(lat_0_rad) / rho)
                 lat_deg = numpy.rad2deg(lat_rad)
                 # fix cases where the distance from the center of the projection is zero
                 lat_deg[rho == 0] = lat_0
                 x1 = x * numpy.sin(c)
                 x2 = rho*numpy.cos(lat_0_rad)*numpy.cos(c) - y*numpy.sin(lat_0_rad)*numpy.sin(c)
                 lon_rad = lon_0_rad + numpy.arctan2(x1, x2)
                 lon_deg = numpy.rad2deg(lon_rad)
                 # Longitudes should be from -180 to 180 degrees
                 lon_deg[lon_deg > 180] -= 360.
                 lon_deg[lon_deg < -180] += 360.
                 return lon_deg, lat_deg
             def antenna_to_geographic(ranges, azimuths, elevations, site):
                 x, y, z = antenna_to_cartesian(numpy.array(ranges), numpy.array(azimuths), numpy.array(elevations))
                 lon, lat = cartesian_to_geographic_aeqd(x, y, site[0], site[1], R=6370997.)
                 return lon, lat
             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 DobleGaussianPlot(SpectraPlot):
                 '''
                 Plot for Double Gaussian Plot
                 '''
                 CODE = 'gaussian_fit'
                 # colormap = 'jet'
                 # plot_type = 'pcolor'
             class DoubleGaussianSpectraCutPlot(SpectraCutPlot):
                 '''
                 Plot SpectraCut with Double Gaussian Fit
                 '''
                 CODE = 'cut_gaussian_fit'
             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/dataOut.normFactor)
                     }
                     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 = 'Range [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 = 10*numpy.log10(self.z)
                     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]
             class WeatherParamsPlot(Plot):
                 #CODE = 'RHI'
                 #plot_name = 'RHI'
                 plot_type = 'scattermap'
                 buffering = False
                 projection = ccrs.PlateCarree()
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.nplots= 1
                     self.ylabel= 'Height [km]'
                     self.xlabel= 'Distance from radar [km]'
                     if self.channels is not None:
                         self.nplots = len(self.channels)
                         self.ncols = len(self.channels)
                     else:
                         self.nplots = self.data.shape(self.CODE)[0]
                         self.ncols = self.nplots
                         self.channels = list(range(self.nplots))
                     self.colorbar=True
                     if len(self.channels)>1:
                         self.width = 12
                     else:
                         self.width   =8
                     self.height  =7
                     self.ini     =0
                     self.len_azi =0
                     self.buffer_ini  = None
                     self.buffer_ele   = None
                     self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.1})
                     self.flag    =0
                     self.indicador= 0
                     self.last_data_ele = None
                     self.val_mean      = None
                 def update(self, dataOut):
                     vars = {
                         'S' : 0,
                         'V' : 1,
                         'W' : 2,
                         'SNR' : 3,
                         'Z' : 4,
                         'D' : 5,
                         'P' : 6,
                         'R' : 7,
                     }
                     data = {}
                     meta = {}
                     if hasattr(dataOut, 'nFFTPoints'):
                         factor = dataOut.normFactor
                     else:
                         factor = 1
                     if hasattr(dataOut, 'dparam'):
                         tmp = getattr(dataOut, 'data_param')
                     else:
+                        #print("-------------------self.attr_data[0]",self.attr_data[0])
                         if 'S' in self.attr_data[0]:
-                            tmp = 10*numpy.log10(10.0*getattr(dataOut, 'data_param')[:,0,:]/(factor))
+                            if self.attr_data[0]=='S':
+                                tmp = 10*numpy.log10(10.0*getattr(dataOut, 'data_param')[:,0,:]/(factor))
+                            if self.attr_data[0]=='SNR':
+                                tmp = 10*numpy.log10(getattr(dataOut, 'data_param')[:,3,:])
                         else:
                             tmp = getattr(dataOut, 'data_param')[:,vars[self.attr_data[0]],:]
                     if self.mask:
                         mask = dataOut.data_param[:,3,:] < self.mask
                         tmp = numpy.ma.masked_array(tmp, mask=mask)
                     r = dataOut.heightList
                     delta_height = r[1]-r[0]
                     valid = numpy.where(r>=0)[0]
                     data['r'] = numpy.arange(len(valid))*delta_height
                     data['data'] = [0, 0]
                     try:
                         data['data'][0] = tmp[0][:,valid]
                         data['data'][1] = tmp[1][:,valid]
                     except:
                         data['data'][0] = tmp[0][:,valid]
                         data['data'][1] = tmp[0][:,valid]
                     if dataOut.mode_op == 'PPI':
                         self.CODE = 'PPI'
                         self.title = self.CODE
                     elif dataOut.mode_op == 'RHI':
                         self.CODE = 'RHI'
                         self.title = self.CODE
                     data['azi'] = dataOut.data_azi
                     data['ele'] = dataOut.data_ele
+                    if isinstance(dataOut.mode_op, bytes):
+                        try:
+                            dataOut.mode_op = dataOut.mode_op.decode()
+                        except:
+                            dataOut.mode_op = str(dataOut.mode_op, 'utf-8')
                     data['mode_op'] = dataOut.mode_op
                     self.mode = dataOut.mode_op
                     return data, meta
                 def plot(self):
                     data = self.data[-1]
                     z = data['data']
                     r = data['r']
                     self.titles = []
                     self.ymax = self.ymax if self.ymax else numpy.nanmax(r)
                     self.ymin = self.ymin if self.ymin else numpy.nanmin(r)
                     self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
                     self.zmin = self.zmin if self.zmin is not None else numpy.nanmin(z)
                     if isinstance(data['mode_op'], bytes):
                         data['mode_op'] = data['mode_op'].decode()
                     if data['mode_op'] == 'RHI':
                         r, theta = numpy.meshgrid(r, numpy.radians(data['ele']))
                         len_aux = int(data['azi'].shape[0]/4)
                         mean = numpy.mean(data['azi'][len_aux:-len_aux])
                         x, y = r*numpy.cos(theta), r*numpy.sin(theta)
                     elif data['mode_op'] == 'PPI':
                         r, theta = numpy.meshgrid(r, -numpy.radians(data['azi'])+numpy.pi/2)
                         len_aux = int(data['ele'].shape[0]/4)
                         mean = numpy.mean(data['ele'][len_aux:-len_aux])
                         x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(mean)), r*numpy.sin(
                                 theta)*numpy.cos(numpy.radians(mean))
                         x = km2deg(x) + -75.295893
                         y = km2deg(y) + -12.040436
                     self.clear_figures()
                     if data['mode_op'] == 'PPI':
                         axes = self.axes['PPI']
                     else:
                         axes = self.axes['RHI']
                     if self.colormap in cb_tables:
                         norm = cb_tables[self.colormap]['norm']
                     else:
                         norm = None
                     for i, ax in enumerate(axes):
                         if data['mode_op'] == 'PPI':
                             ax.set_extent([-75.745893, -74.845893, -12.490436, -11.590436])
                         if norm is None:
                             ax.plt = ax.pcolormesh(x, y, z[i], cmap=self.colormap, vmin=self.zmin, vmax=self.zmax)
                         else:
                             ax.plt = ax.pcolormesh(x, y, z[i], cmap=self.colormap, norm=norm)
                         if data['mode_op'] == 'RHI':
                             len_aux = int(data['azi'].shape[0]/4)
                             mean = numpy.mean(data['azi'][len_aux:-len_aux])
                             if len(self.channels) !=1:
                                 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels]
                             else:
                                 self.titles = ['RHI {} at AZ: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
                         elif data['mode_op'] == 'PPI':
                             len_aux = int(data['ele'].shape[0]/4)
                             mean = numpy.mean(data['ele'][len_aux:-len_aux])
                             if len(self.channels) !=1:
                                 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[x], str(round(mean,1)), x) for x in self.channels]
                             else:
                                 self.titles = ['PPI {} at EL: {} CH {}'.format(self.labels[0], str(round(mean,1)), self.channels[0])]
                         self.mode_value = round(mean,1)
                         if data['mode_op'] == 'PPI':
                             gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                               linewidth=1, color='gray', alpha=0.5, linestyle='--')
                             gl.xlabel_style = {'size': 8}
                             gl.ylabel_style = {'size': 8}
                             gl.xlabels_top = False
                             gl.ylabels_right = False
+                            #self.shapes="/home/soporte/workspace/sirm/volumes/schain/shapes/"
+                            #print("self.shapes",self.shapes)
                             shape_p = os.path.join(self.shapes,'PER_ADM2/PER_ADM2.shp')
                             shape_d = os.path.join(self.shapes,'PER_ADM1/PER_ADM1.shp')
                             capitales = os.path.join(self.shapes,'CAPITALES/cap_provincia.shp')
                             vias = os.path.join(self.shapes,'Carreteras/VIAS_NACIONAL_250000.shp')
                             reader_d = shpreader.BasicReader(shape_p, encoding='latin1')
                             reader_p = shpreader.BasicReader(shape_d, encoding='latin1')
                             reader_c = shpreader.BasicReader(capitales, encoding='latin1')
                             reader_v = shpreader.BasicReader(vias, encoding='latin1')
                             caps = [x for x in reader_c.records()  if x.attributes["Departa"] in ("JUNIN", "LIMA", "AYACUCHO", "HUANCAVELICA")]
                             districts = [x for x in reader_d.records() if x.attributes["Name"] in ("JUNÍN", "CHANCHAMAYO", "CHUPACA", "CONCEPCIÓN", "HUANCAYO", "JAUJA", "SATIPO", "TARMA", "YAUYOS", "HUAROCHIRÍ", "CANTA", "HUANTA", "TAYACAJA")]
                             provs = [x for x in reader_p.records() if x.attributes["NAME"] in ("Junín", "Lima")]
                             vias = [x for x in reader_v.records() if x.attributes["DEP"] in ("JUNIN", "LIMA")]
                             # Display limits and streets
                             shape_feature = ShapelyFeature([x.geometry for x in districts], ccrs.PlateCarree(), facecolor="none", edgecolor='grey', lw=0.5)
                             ax.add_feature(shape_feature)
                             shape_feature = ShapelyFeature([x.geometry for x in provs], ccrs.PlateCarree(), facecolor="none", edgecolor='white', lw=1)
                             ax.add_feature(shape_feature)
                             shape_feature = ShapelyFeature([x.geometry for x in vias], ccrs.PlateCarree(), facecolor="none", edgecolor='yellow', lw=1)
                             ax.add_feature(shape_feature)
                             for cap in caps:
                                 if cap.attributes['Nombre'] in ("LA OROYA", "CONCEPCIÓN", "HUANCAYO", "JAUJA", "CHUPACA", "YAUYOS", "HUANTA", "PAMPAS"):
                                     ax.text(cap.attributes['X'], cap.attributes['Y'], cap.attributes['Nombre'].title(), size=7, color='white')
                             ax.text(-75.052003, -11.915552, 'Huaytapallana', size=7, color='cyan')
                             ax.plot(-75.052003, -11.915552, '*')
                             for R in (10, 20, 30 , 40, 50):
                                 circle = Circle((-75.295893, -12.040436), km2deg(R), facecolor='none',
                                     edgecolor='skyblue', linewidth=1, alpha=0.5)
                                 ax.add_patch(circle)
                                 ax.text(km2deg(R)*numpy.cos(numpy.radians(45))-75.295893,
                                     km2deg(R)*numpy.sin(numpy.radians(45))-12.040436,
                                     '{}km'.format(R), color='skyblue', size=7)
                         elif data['mode_op'] == 'RHI':
                             ax.grid(color='grey', alpha=0.5, linestyle='--', linewidth=1)

schainpy/model/graphics/jroplot_spectra.py +1 -1

             # Copyright (c) 2012-2021 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
                 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):
                     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))
+                    meta['xrange'] = (dataOut.getFreqRange(0)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(0))
                     if self.CODE == 'spc_moments':
                         data['moments'] = dataOut.moments
                         # data['spc'] = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
                     if self.CODE == 'gaussian_fit':
                         # data['moments'] = dataOut.moments
                         data['gaussfit'] = dataOut.DGauFitParams
                         # data['spc'] = 10*numpy.log10(dataOut.data_pre[0]/dataOut.normFactor)
                     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') | (self.CODE == 'gaussian_fit'):
                         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 self.CODE == 'gaussian_fit':
                             # mean = data['moments'][n, 1]
                             gau0 = data['gaussfit'][n][2,:,0]
                             gau1 = data['gaussfit'][n][2,:,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', lw=1)[0]
                             if self.CODE == 'gaussian_fit':
                                 # ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
                                 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
                                 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[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)
                             if self.CODE == 'gaussian_fit':
                                 # ax.plt_mean.set_data(mean, y)
                                 ax.plt_gau0.set_data(gau0, y)
                                 ax.plt_gau1.set_data(gau1, y)
                         self.titles.append('CH {}: {:3.2f}dB'.format(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'
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     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.nrows)]
                 def update(self, dataOut):
                     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.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_pre[0]/dataOut.normFactor)
                     data['spc'] = spc
                     meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
                     if self.CODE == 'cut_gaussian_fit':
                         data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
                         data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
                     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][:-1]
                         self.xlabel = "Velocity (m/s)"
                     if self.CODE == 'cut_gaussian_fit':
                         x = self.data.xrange[2][:-1]
                         self.xlabel = "Velocity (m/s)"
                     self.titles = []
                     y = self.data.yrange
                     data = self.data[-1]
                     z = data['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 self.CODE == 'cut_gaussian_fit':
                             gau0 = data['gauss_fit0']
                             gau1 = data['gauss_fit1']
                         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, lw=0.25)
                             if self.CODE == 'cut_gaussian_fit':
                                 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
                                 for i, line in enumerate(ax.plt_gau0):
                                     line.set_color(ax.plt[i].get_color())
                                 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
                                 for i, line in enumerate(ax.plt_gau1):
                                     line.set_color(ax.plt[i].get_color())
                             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]].T)
                             for i, line in enumerate(ax.plt_gau0):
                                 line.set_data(x, gau0[n, :, index[i]].T)
                                 line.set_color(ax.plt[i].get_color())
                             for i, line in enumerate(ax.plt_gau1):
                                 line.set_data(x, gau1[n, :, index[i]].T)
                                 line.set_color(ax.plt[i].get_color())
                         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_digitalRF.py +33 -14

             '''
             Created on Jul 3, 2014
             @author: roj-idl71
             '''
             # SUBCHANNELS EN VEZ DE CHANNELS
             # BENCHMARKS -> PROBLEMAS CON ARCHIVOS GRANDES -> INCONSTANTE EN EL TIEMPO
             # ACTUALIZACION DE VERSION
             # HEADERS
             # MODULO DE ESCRITURA
             # METADATA
             import os
             import time
             import datetime
             import numpy
             import timeit
             from fractions import Fraction
             from time import time
             from time import sleep
             import schainpy.admin
             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
             import pickle
             try:
+                os.environ["HDF5_USE_FILE_LOCKING"] = "FALSE"
                 import digital_rf
             except:
                 pass
             class DigitalRFReader(ProcessingUnit):
                 '''
                 classdocs
                 '''
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     ProcessingUnit.__init__(self)
                     self.dataOut                  = Voltage()
                     self.__printInfo              = True
                     self.__flagDiscontinuousBlock = False
                     self.__bufferIndex = 9999999
                     self.__codeType    = 0
                     self.__ippKm       = None
                     self.__nCode       = None
                     self.__nBaud       = None
                     self.__code        = None
                     self.dtype         = None
                     self.oldAverage    = None
                     self.path          = None
                 def close(self):
                     print('Average of writing to digital rf format is ', self.oldAverage * 1000)
                     return
                 def __getCurrentSecond(self):
                     return self.__thisUnixSample / self.__sample_rate
                 thisSecond = property(__getCurrentSecond, "I'm the 'thisSecond' property.")
                 def __setFileHeader(self):
                     '''
                     In this method will be initialized every parameter of dataOut object (header, no data)
                     '''
                     ippSeconds = 1.0 * self.__nSamples / self.__sample_rate
                     if not self.getByBlock:
                         nProfiles = 1.0 / ippSeconds  # Number of profiles in one second
                     else:
                         nProfiles = self.nProfileBlocks # Number of profiles in one block
                     try:
                         self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
                             self.__radarControllerHeader)
                     except:
                         self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
                             txA=0,
                             txB=0,
                             nWindows=1,
                             nHeights=self.__nSamples,
                             firstHeight=self.__firstHeigth,
                             deltaHeight=self.__deltaHeigth,
                             codeType=self.__codeType,
                             nCode=self.__nCode, nBaud=self.__nBaud,
                             code=self.__code)
                     try:
                         self.dataOut.systemHeaderObj = SystemHeader(self.__systemHeader)
                     except:
                         self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
                                                                     nProfiles=nProfiles,
                                                                     nChannels=len(
                                                                         self.__channelList),
                                                                     adcResolution=14)
                     self.dataOut.type  = "Voltage"
                     self.dataOut.data  = None
                     self.dataOut.dtype = self.dtype
                     # self.dataOut.nChannels = 0
                     # self.dataOut.nHeights = 0
                     self.dataOut.nProfiles   = int(nProfiles)
                     self.dataOut.heightList  = self.__firstHeigth + \
                         numpy.arange(self.__nSamples, dtype=numpy.float) * \
                         self.__deltaHeigth
                     #self.dataOut.channelList = list(range(self.__num_subchannels))
                     self.dataOut.channelList = list(range(len(self.__channelList)))
                     if not self.getByBlock:
                         self.dataOut.blocksize   = self.dataOut.nChannels * self.dataOut.nHeights
                     else:
                         self.dataOut.blocksize   = self.dataOut.nChannels * self.dataOut.nHeights*self.nProfileBlocks
                     # self.dataOut.channelIndexList = None
                     self.dataOut.flagNoData      = True
                     if not self.getByBlock:
                         self.dataOut.flagDataAsBlock = False
                     else:
                         self.dataOut.flagDataAsBlock = True
                     # Set to TRUE if the data is discontinuous
                     self.dataOut.flagDiscontinuousBlock = False
                     self.dataOut.utctime     = None
                     # timezone like jroheader, difference in minutes between UTC and localtime
                     self.dataOut.timeZone    = self.__timezone / 60
                     self.dataOut.dstFlag     = 0
                     self.dataOut.errorCount  = 0
                     try:
                         self.dataOut.nCohInt = self.fixed_metadata_dict.get(
                             'nCohInt', self.nCohInt)
                         # asumo que la data esta decodificada
                         self.dataOut.flagDecodeData = self.fixed_metadata_dict.get(
                             'flagDecodeData', self.flagDecodeData)
                         # asumo que la data esta sin flip
                         self.dataOut.flagDeflipData = self.fixed_metadata_dict['flagDeflipData']
                         self.dataOut.flagShiftFFT   = self.fixed_metadata_dict['flagShiftFFT']
                         self.dataOut.useLocalTime   = self.fixed_metadata_dict['useLocalTime']
                     except:
                         pass
                     self.dataOut.ippSeconds = ippSeconds
                     # Time interval between profiles
                     # self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
                     self.dataOut.frequency  = self.__frequency
                     self.dataOut.realtime   = self.__online
                 def findDatafiles(self, path, startDate=None, endDate=None):
                     if not os.path.isdir(path):
                         return []
                     try:
                         digitalReadObj = digital_rf.DigitalRFReader(
                             path, load_all_metadata=True)
                     except:
                         digitalReadObj = digital_rf.DigitalRFReader(path)
                     channelNameList    = digitalReadObj.get_channels()
                     if not channelNameList:
                         return []
                     metadata_dict      = digitalReadObj.get_rf_file_metadata(channelNameList[0])
                     sample_rate        = metadata_dict['sample_rate'][0]
                     this_metadata_file = digitalReadObj.get_metadata(channelNameList[0])
                     try:
                         timezone       = this_metadata_file['timezone'].value
                     except:
                         timezone       = 0
                     startUTCSecond, endUTCSecond = digitalReadObj.get_bounds(
                         channelNameList[0]) / sample_rate - timezone
                     startDatetime      = datetime.datetime.utcfromtimestamp(startUTCSecond)
                     endDatatime        = datetime.datetime.utcfromtimestamp(endUTCSecond)
                     if not startDate:
                         startDate      = startDatetime.date()
                     if not endDate:
                         endDate        = endDatatime.date()
                     dateList           = []
                     thisDatetime       = startDatetime
                     while(thisDatetime <= endDatatime):
                         thisDate        = thisDatetime.date()
                         if thisDate     < startDate:
                             continue
                         if thisDate     > endDate:
                             break
                         dateList.append(thisDate)
                         thisDatetime += datetime.timedelta(1)
                     return dateList
                 def setup(self, path=None,
                           startDate=None,
                           endDate=None,
                           startTime=datetime.time(0, 0, 0),
                           endTime=datetime.time(23, 59, 59),
                           channelList=None,
                           nSamples=None,
                           online=False,
                           delay=60,
                           buffer_size=1024,
                           ippKm=None,
                           nCohInt=1,
                           nCode=1,
                           nBaud=1,
                           flagDecodeData=False,
                           code=numpy.ones((1, 1), dtype=numpy.int),
                           getByBlock=0,
                           nProfileBlocks=1,
                           **kwargs):
                     '''
                     In this method we should set all initial parameters.
                     Inputs:
                         path
                         startDate
                         endDate
                         startTime
                         endTime
                         set
                         expLabel
                         ext
                         online
                         delay
                     '''
                     self.path           = path
                     self.nCohInt        = nCohInt
                     self.flagDecodeData = flagDecodeData
                     self.i              = 0
                     self.getByBlock     = getByBlock
                     self.nProfileBlocks = nProfileBlocks
+                    if online:
+                        print('Waiting for RF data..')
+                        sleep(40)
                     if not os.path.isdir(path):
                         raise ValueError("[Reading] Directory %s does not exist" % path)
+                    #print("path",path)
                     try:
                         self.digitalReadObj = digital_rf.DigitalRFReader(
                             path, load_all_metadata=True)
                     except:
                         self.digitalReadObj = digital_rf.DigitalRFReader(path)
                     channelNameList         = self.digitalReadObj.get_channels()
                     if not channelNameList:
                         raise ValueError("[Reading] Directory %s does not have any files" % path)
                     if not channelList:
                         channelList = list(range(len(channelNameList)))
                     ##########  Reading metadata ######################
                     top_properties           = self.digitalReadObj.get_properties(
                         channelNameList[channelList[0]])
                     self.__num_subchannels   = top_properties['num_subchannels']
                     self.__sample_rate       = 1.0 * \
                         top_properties['sample_rate_numerator'] / \
                         top_properties['sample_rate_denominator']
                     # self.__samples_per_file = top_properties['samples_per_file'][0]
                     self.__deltaHeigth       = 1e6 * 0.15 / self.__sample_rate  # why 0.15?
                     this_metadata_file       = self.digitalReadObj.get_digital_metadata(
                         channelNameList[channelList[0]])
                     metadata_bounds          = this_metadata_file.get_bounds()
                     self.fixed_metadata_dict = this_metadata_file.read(
                         metadata_bounds[0])[metadata_bounds[0]]  # GET FIRST HEADER
                     try:
                         self.__processingHeader      = self.fixed_metadata_dict['processingHeader']
                         self.__radarControllerHeader = self.fixed_metadata_dict['radarControllerHeader']
                         self.__systemHeader          = self.fixed_metadata_dict['systemHeader']
                         self.dtype                   = pickle.loads(self.fixed_metadata_dict['dtype'])
                     except:
                         pass
                     self.__frequency = None
                     self.__frequency = self.fixed_metadata_dict.get('frequency', 1)
                     self.__timezone = self.fixed_metadata_dict.get('timezone', 18000)
                     try:
                         nSamples = self.fixed_metadata_dict['nSamples']
                     except:
                         nSamples = None
                     self.__firstHeigth = 0
                     try:
                         codeType       = self.__radarControllerHeader['codeType']
                     except:
                         codeType       = 0
                     try:
                         if codeType:
                             nCode = self.__radarControllerHeader['nCode']
                             nBaud = self.__radarControllerHeader['nBaud']
                             code  = self.__radarControllerHeader['code']
                     except:
                         pass
                     if not ippKm:
                         try:
                             # seconds to km
                             ippKm = self.__radarControllerHeader['ipp']
                         except:
                             ippKm = None
                     ####################################################
                     self.__ippKm   = ippKm
                     startUTCSecond = None
                     endUTCSecond   = None
                     if startDate:
                         startDatetime  = datetime.datetime.combine(startDate, startTime)
                         startUTCSecond = (
-                            startDatetime - datetime.datetime(1970, 1, 1)).total_seconds() + self.__timezone
+                            startDatetime - datetime.datetime(1970, 1, 1)).total_seconds()# + self.__timezone
                     if endDate:
                         endDatetime   = datetime.datetime.combine(endDate, endTime)
                         endUTCSecond  = (endDatetime - datetime.datetime(1970,
-, 1)).total_seconds() + self.__timezone
+, 1)).total_seconds()# + self.__timezone
+                    start_index, end_index = self.digitalReadObj.get_bounds(channelNameList[channelList[0]])
+                    if start_index==None or end_index==None:
-                    #print(startUTCSecond,endUTCSecond)
+                         print("Check error No data,  start_index: ",start_index,",end_index: ",end_index)
-                    start_index, end_index = self.digitalReadObj.get_bounds(
+                         #return 0
-                        channelNameList[channelList[0]])
-                    #print("*****",start_index,end_index)
                     if not startUTCSecond:
                         startUTCSecond = start_index / self.__sample_rate
                     if start_index     > startUTCSecond * self.__sample_rate:
                         startUTCSecond = start_index / self.__sample_rate
                     if not endUTCSecond:
                         endUTCSecond   = end_index / self.__sample_rate
                     if end_index       < endUTCSecond * self.__sample_rate:
                         endUTCSecond   = end_index / self.__sample_rate #Check UTC and LT time
                     if not nSamples:
                         if not ippKm:
                             raise ValueError("[Reading] nSamples or ippKm should be defined")
                         nSamples            = int(ippKm / (1e6 * 0.15 / self.__sample_rate))
                     channelBoundList        = []
                     channelNameListFiltered = []
                     for thisIndexChannel in channelList:
                         thisChannelName        = channelNameList[thisIndexChannel]
                         start_index, end_index = self.digitalReadObj.get_bounds(
                             thisChannelName)
                         channelBoundList.append((start_index, end_index))
                         channelNameListFiltered.append(thisChannelName)
                     self.profileIndex = 0
                     self.i            = 0
                     self.__delay      = delay
                     self.__codeType   = codeType
                     self.__nCode      = nCode
                     self.__nBaud      = nBaud
                     self.__code       = code
                     self.__datapath         = path
                     self.__online           = online
                     self.__channelList      = channelList
                     self.__channelNameList  = channelNameListFiltered
                     self.__channelBoundList = channelBoundList
                     self.__nSamples         = nSamples
                     if self.getByBlock:
                         nSamples = nSamples*nProfileBlocks
                     self.__samples_to_read  = int(nSamples)  # FIJO: AHORA 40
                     self.__nChannels        = len(self.__channelList)
                     #print("------------------------------------------")
                     #print("self.__samples_to_read",self.__samples_to_read)
                     #print("self.__nSamples",self.__nSamples)
                     # son iguales y el buffer_index da 0
                     self.__startUTCSecond   = startUTCSecond
                     self.__endUTCSecond     = endUTCSecond
                     self.__timeInterval     = 1.0 * self.__samples_to_read / \
                         self.__sample_rate  # Time interval
                     if online:
                         # self.__thisUnixSample = int(endUTCSecond*self.__sample_rate - 4*self.__samples_to_read)
                         startUTCSecond = numpy.floor(endUTCSecond)
                     # por que en el otro metodo lo primero q se hace es sumar samplestoread
                     self.__thisUnixSample = int(startUTCSecond * self.__sample_rate) - self.__samples_to_read
                     #self.__data_buffer    = numpy.zeros(
                     #    (self.__num_subchannels, self.__samples_to_read), dtype=numpy.complex)
+                    print("samplestoread",self.__samples_to_read)
                     self.__data_buffer    = numpy.zeros((int(len(channelList)), self.__samples_to_read), dtype=numpy.complex)
                     self.__setFileHeader()
                     self.isConfig = True
                     print("[Reading] Digital RF Data was found from %s to %s " % (
                         datetime.datetime.utcfromtimestamp(
                             self.__startUTCSecond - self.__timezone),
                         datetime.datetime.utcfromtimestamp(
                             self.__endUTCSecond - self.__timezone)
                     ))
                     print("[Reading] Starting process from %s to %s" % (datetime.datetime.utcfromtimestamp(startUTCSecond - self.__timezone),
-                                                                        datetime.datetime.utcfromtimestamp(
+                                                                        datetime.datetime.utcfromtimestamp(endUTCSecond - self.__timezone)))
-                        endUTCSecond - self.__timezone)
-                    ))
                     self.oldAverage    = None
                     self.count         = 0
                     self.executionTime = 0
                 def __reload(self):
                     #         print
                     #         print "%s not in range [%s, %s]" %(
                     #                                           datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
                     #                                           datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone),
                     #                                           datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone)
                     #                                           )
                     print("[Reading] reloading metadata ...")
                     try:
                         self.digitalReadObj.reload(complete_update=True)
                     except:
                         self.digitalReadObj = digital_rf.DigitalRFReader(self.path)
                     start_index, end_index  = self.digitalReadObj.get_bounds(
                         self.__channelNameList[self.__channelList[0]])
                     if start_index          > self.__startUTCSecond * self.__sample_rate:
                         self.__startUTCSecond = 1.0 * start_index / self.__sample_rate
                     if end_index            > self.__endUTCSecond * self.__sample_rate:
                         self.__endUTCSecond = 1.0 * end_index / self.__sample_rate
                         print()
                         print("[Reading] New timerange found [%s, %s] " % (
                             datetime.datetime.utcfromtimestamp(
                                 self.__startUTCSecond - self.__timezone),
                             datetime.datetime.utcfromtimestamp(
                                 self.__endUTCSecond - self.__timezone)
                         ))
                         return True
                     return False
                 def timeit(self, toExecute):
                     t0                  = time.time()
                     toExecute()
                     self.executionTime  = time.time() - t0
                     if self.oldAverage is None:
                         self.oldAverage = self.executionTime
                     self.oldAverage     = (self.executionTime + self.count *
                                        self.oldAverage) / (self.count + 1.0)
                     self.count          = self.count + 1.0
                     return
                 def __readNextBlock(self, seconds=30, volt_scale=1/20000.0):
                     '''
                     NOTA: APLICACION RADAR METEOROLOGICO
                     VALORES OBTENIDOS CON LA USRP, volt_scale = 1,conexion directa al Ch Rx.
                     MAXIMO
 -> 0.980 Voltiospp
 -> 0.480 Voltiospp
 -> 1.440 Voltiospp
 -> 1.940 Voltiospp
                     Para llevar al valor correspondiente de Voltaje, debemos dividir por 20000
                     y obtenemos la Amplitud correspondiente de entrada IQ.
                     volt_scale = (1/20000.0)
                     '''
                     # Set the next data
                     self.__flagDiscontinuousBlock = False
                     self.__thisUnixSample        += self.__samples_to_read
                     if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
                         print ("[Reading] There are no more data into selected time-range")
                         if self.__online:
                             sleep(3)
                             self.__reload()
                         else:
                             return False
                         if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
                             return False
                             self.__thisUnixSample -= self.__samples_to_read
                     indexChannel = 0
                     dataOk = False
                     for thisChannelName in self.__channelNameList:  # TODO VARIOS CHANNELS?
                         for indexSubchannel in range(self.__num_subchannels):
                             try:
                                 t0     = time()
+                                #print("thisUNixSample",self.__thisUnixSample)
                                 result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample,
                                                                               self.__samples_to_read,
                                                                               thisChannelName, sub_channel=indexSubchannel)
+                                #print("result--------------",result)
                                 self.executionTime  = time() - t0
                                 if self.oldAverage is None:
                                     self.oldAverage = self.executionTime
                                 self.oldAverage     = (
                                     self.executionTime + self.count * self.oldAverage) / (self.count + 1.0)
                                 self.count = self.count + 1.0
                             except IOError as e:
                                 # read next profile
                                 self.__flagDiscontinuousBlock = True
                                 print("[Reading] %s" % datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), e)
-                                break
+                                bot = 0
+                                while(self.__flagDiscontinuousBlock):
+                                    bot +=1
+                                    self.__thisUnixSample += self.__sample_rate
+                                    try:
+                                        result = result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample,self.__samples_to_read,thisChannelName, sub_channel=indexSubchannel)
+                                        self.__flagDiscontinuousBlock=False
+                                        print("Searching.. N°: ",bot,"Success",self.__thisUnixSample)
+                                    except:
+                                        print("Searching...N°: ",bot,"Fail", self.__thisUnixSample)
+                                if self.__flagDiscontinuousBlock==True:
+                                    break
+                                else:
+                                    print("New data index found...",self.__thisUnixSample)
+                                #break
                             if result.shape[0] != self.__samples_to_read:
                                 self.__flagDiscontinuousBlock = True
                                 print("[Reading] %s: Too few samples were found, just %d/%d  samples" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
                                                                                                          result.shape[0],
                                                                                                          self.__samples_to_read))
                                 break
                             self.__data_buffer[indexChannel, :] = result * volt_scale
                             indexChannel+=1
                             dataOk       = True
                     self.__utctime       = self.__thisUnixSample / self.__sample_rate
                     if not dataOk:
                         return False
                     print("[Reading] %s: %d samples <> %f sec" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
                                                                   self.__samples_to_read,
                                                                   self.__timeInterval))
                     self.__bufferIndex   = 0
                     return True
                 def __isBufferEmpty(self):
                     return self.__bufferIndex > self.__samples_to_read - self.__nSamples  # 40960 - 40
                 def getData(self, seconds=30, nTries=5):
                     '''
                         This method gets the data from files and put the data into the dataOut object
                         In addition, increase el the buffer counter in one.
                         Return:
                             data    :    retorna un perfil de voltages (alturas * canales) copiados desde el
                                         buffer. Si no hay mas archivos a leer retorna None.
                         Affected:
                             self.dataOut
                             self.profileIndex
                             self.flagDiscontinuousBlock
                             self.flagIsNewBlock
                     '''
                     #print("getdata")
                     err_counter = 0
                     self.dataOut.flagNoData = True
                     if self.__isBufferEmpty():
                         #print("hi")
                         self.__flagDiscontinuousBlock = False
                         while True:
                             if self.__readNextBlock():
                                 break
                             if self.__thisUnixSample > self.__endUTCSecond * self.__sample_rate:
                                 raise schainpy.admin.SchainError('Error')
                                 return
                             if self.__flagDiscontinuousBlock:
                                 raise schainpy.admin.SchainError('discontinuous block found')
                                 return
                             if not self.__online:
                                 raise schainpy.admin.SchainError('Online?')
                                 return
                             err_counter += 1
                             if err_counter > nTries:
                                 raise schainpy.admin.SchainError('Max retrys reach')
                                 return
                             print('[Reading] waiting %d seconds to read a new block' % seconds)
                             sleep(seconds)
                         if not self.getByBlock:
                             #print("self.__bufferIndex",self.__bufferIndex)# este valor siempre es cero aparentemente
                             self.dataOut.data = self.__data_buffer[:, self.__bufferIndex:self.__bufferIndex + self.__nSamples]
                             self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
                             self.dataOut.flagNoData = False
                             self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
                             self.dataOut.profileIndex = self.profileIndex
                             self.__bufferIndex += self.__nSamples
                             self.profileIndex  += 1
                             if self.profileIndex == self.dataOut.nProfiles:
                                 self.profileIndex = 0
                         else:
                             # ojo debo anadir el readNextBLock y el  __isBufferEmpty(
                             self.dataOut.flagNoData             = False
                             buffer = self.__data_buffer[:,self.__bufferIndex:self.__bufferIndex + self.__samples_to_read]
                             buffer = buffer.reshape((self.__nChannels, self.nProfileBlocks, int(self.__samples_to_read/self.nProfileBlocks)))
                             self.dataOut.nProfileBlocks = self.nProfileBlocks
                             self.dataOut.data = buffer
                             self.dataOut.utctime = ( self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
                             self.profileIndex  += self.__samples_to_read
                             self.__bufferIndex += self.__samples_to_read
                             self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
                         return True
                 def printInfo(self):
                     '''
                     '''
                     if self.__printInfo == False:
                         return
                     # self.systemHeaderObj.printInfo()
                     # self.radarControllerHeaderObj.printInfo()
                     self.__printInfo = False
                 def printNumberOfBlock(self):
                     '''
                     '''
                     return
                     # print self.profileIndex
                 def run(self, **kwargs):
                     '''
                     This method will be called many times so here you should put all your code
                     '''
                     if not self.isConfig:
                         self.setup(**kwargs)
                     self.getData(seconds=self.__delay)
                     return
             @MPDecorator
             class DigitalRFWriter(Operation):
                 '''
                 classdocs
                 '''
                 def __init__(self, **kwargs):
                     '''
                     Constructor
                     '''
                     Operation.__init__(self, **kwargs)
                     self.metadata_dict = {}
                     self.dataOut       = None
                     self.dtype         = None
                     self.oldAverage    = 0
                 def setHeader(self):
                     self.metadata_dict['frequency']      = self.dataOut.frequency
                     self.metadata_dict['timezone']       = self.dataOut.timeZone
                     self.metadata_dict['dtype']          = pickle.dumps(self.dataOut.dtype)
                     self.metadata_dict['nProfiles']      = self.dataOut.nProfiles
                     self.metadata_dict['heightList']     = self.dataOut.heightList
                     self.metadata_dict['channelList']    = self.dataOut.channelList
                     self.metadata_dict['flagDecodeData'] = self.dataOut.flagDecodeData
                     self.metadata_dict['flagDeflipData'] = self.dataOut.flagDeflipData
                     self.metadata_dict['flagShiftFFT']   = self.dataOut.flagShiftFFT
                     self.metadata_dict['useLocalTime']   = self.dataOut.useLocalTime
                     self.metadata_dict['nCohInt']        = self.dataOut.nCohInt
                     self.metadata_dict['type']           = self.dataOut.type
                     self.metadata_dict['flagDataAsBlock']= getattr(
                         self.dataOut, 'flagDataAsBlock', None)  # chequear
                 def setup(self, dataOut, path, frequency, fileCadence, dirCadence, metadataCadence, set=0, metadataFile='metadata', ext='.h5'):
                     '''
                     In this method we should set all initial parameters.
                     Input:
                         dataOut: Input data will also be outputa data
                     '''
                     self.setHeader()
                     self.__ippSeconds  = dataOut.ippSeconds
                     self.__deltaH      = dataOut.getDeltaH()
                     self.__sample_rate = 1e6 * 0.15 / self.__deltaH
                     self.__dtype       = dataOut.dtype
                     if len(dataOut.dtype) == 2:
                         self.__dtype = dataOut.dtype[0]
                     self.__nSamples  = dataOut.systemHeaderObj.nSamples
                     self.__nProfiles = dataOut.nProfiles
                     if self.dataOut.type != 'Voltage':
                         raise 'Digital RF cannot be used with this data type'
                         self.arr_data = numpy.ones((1, dataOut.nFFTPoints * len(
                             self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
                     else:
                         self.arr_data = numpy.ones((self.__nSamples, len(
                             self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
                     file_cadence_millisecs  = 1000
                     sample_rate_fraction    = Fraction(self.__sample_rate).limit_denominator()
                     sample_rate_numerator   = int(sample_rate_fraction.numerator)
                     sample_rate_denominator = int(sample_rate_fraction.denominator)
                     start_global_index      = dataOut.utctime * self.__sample_rate
                     uuid              = 'prueba'
                     compression_level = 0
                     checksum          = False
                     is_complex        = True
                     num_subchannels   = len(dataOut.channelList)
                     is_continuous     = True
                     marching_periods  = False
                     self.digitalWriteObj = digital_rf.DigitalRFWriter(path, self.__dtype, dirCadence,
                                                                       fileCadence, start_global_index,
                                                                       sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum,
                                                                       is_complex, num_subchannels, is_continuous, marching_periods)
                     metadata_dir         = os.path.join(path, 'metadata')
                     os.system('mkdir %s' % (metadata_dir))
                     self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1,  # 236, file_cadence_millisecs / 1000
                                                                                     sample_rate_numerator, sample_rate_denominator,
                                                                                     metadataFile)
                     self.isConfig      = True
                     self.currentSample = 0
                     self.oldAverage    = 0
                     self.count         = 0
                     return
                 def writeMetadata(self):
                     start_idx                                   = self.__sample_rate * self.dataOut.utctime
                     self.metadata_dict['processingHeader']      = self.dataOut.processingHeaderObj.getAsDict(
                     )
                     self.metadata_dict['radarControllerHeader'] = self.dataOut.radarControllerHeaderObj.getAsDict(
                     )
                     self.metadata_dict['systemHeader']          = self.dataOut.systemHeaderObj.getAsDict(
                     )
                     self.digitalMetadataWriteObj.write(start_idx, self.metadata_dict)
                     return
                 def timeit(self, toExecute):
                     t0 = time()
                     toExecute()
                     self.executionTime  = time() - t0
                     if self.oldAverage is None:
                         self.oldAverage = self.executionTime
                     self.oldAverage     = (self.executionTime + self.count *
                                        self.oldAverage) / (self.count + 1.0)
                     self.count          = self.count + 1.0
                     return
                 def writeData(self):
                     if self.dataOut.type != 'Voltage':
                         raise 'Digital RF cannot be used with this data type'
                         for channel in self.dataOut.channelList:
                             for i in range(self.dataOut.nFFTPoints):
                                 self.arr_data[1][channel * self.dataOut.nFFTPoints +
                                                  i]['r'] = self.dataOut.data[channel][i].real
                                 self.arr_data[1][channel * self.dataOut.nFFTPoints +
                                                  i]['i'] = self.dataOut.data[channel][i].imag
                     else:
                         for i in range(self.dataOut.systemHeaderObj.nSamples):
                             for channel in self.dataOut.channelList:
                                 self.arr_data[i][channel]['r'] = self.dataOut.data[channel][i].real
                                 self.arr_data[i][channel]['i'] = self.dataOut.data[channel][i].imag
                     def f(): return self.digitalWriteObj.rf_write(self.arr_data)
                     self.timeit(f)
                     return
                 def run(self, dataOut, frequency=49.92e6, path=None, fileCadence=1000, dirCadence=36000, metadataCadence=1, **kwargs):
                     '''
                     This method will be called many times so here you should put all your code
                     Inputs:
                         dataOut: object with the data
                     '''
                     # print dataOut.__dict__
                     self.dataOut = dataOut
                     if not self.isConfig:
                         self.setup(dataOut, path, frequency, fileCadence,
                                    dirCadence, metadataCadence, **kwargs)
                         self.writeMetadata()
                     self.writeData()
                     ## self.currentSample += 1
                     # if self.dataOut.flagDataAsBlock or self.currentSample == 1:
                     # self.writeMetadata()
                     ## if self.currentSample == self.__nProfiles: self.currentSample = 0
                     return dataOut# en la version 2.7 no aparece este return
                 def close(self):
                     print('[Writing] - Closing files ')
                     print('Average of writing to digital rf format is ', self.oldAverage * 1000)
                     try:
                         self.digitalWriteObj.close()
                     except:
                         pass

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/scripts/sophy_proc.py +2 -2

             # SOPHY PROC script
             import os, sys, json, argparse
             import datetime
             import time
             PATH = '/DATA_RM/DATA'
             PATH = '/media/jespinoza/Elements'
             PATH = '/media/jespinoza/data/SOPHY'
             PATH = '/home/soporte/Documents/EVENTO'
             PARAM = {
                 'S': {'zmin': -45, 'zmax': -25, 'colormap': 'jet', 'label': 'Power', 'wrname': 'power','cb_label': 'dBm', 'ch':0},
                 'SNR': {'zmin': -40, 'zmax': -20, 'colormap': 'jet', 'label': 'SNR', 'wrname': 'snr','cb_label': 'dB', 'ch':0},
                 'V': {'zmin': -12, 'zmax': 12, 'colormap': 'sophy_v', 'label': 'Velocity', 'wrname': 'velocity', 'cb_label': 'm/s', 'ch':0},
                 'R': {'zmin': 0,   'zmax': 1,  'colormap': 'jet',    'label': 'RhoHV', 'wrname':'rhoHV', 'cb_label': '*',  'ch':0},
                 'P': {'zmin': -180,'zmax': 180,'colormap': 'RdBu_r', 'label': 'PhiDP', 'wrname':'phiDP' , 'cb_label': 'º',  'ch':0},
                 'D': {'zmin': -30, 'zmax': 80, 'colormap': 'sophy_r','label': 'ZDR','wrname':'differential_reflectivity' , 'cb_label': 'dBz','ch':0},
                 'Z':  {'zmin': -30, 'zmax': 80, 'colormap': 'sophy_r','label': 'Reflectivity ',  'wrname':'reflectivity', 'cb_label': 'dBz','ch':0},
                 'W':  {'zmin': 0, 'zmax': 15, 'colormap': 'sophy_w','label': 'Spectral Width', 'wrname':'spectral_width', 'cb_label': 'm/s', 'ch':0}
                 }
             def max_index(r, sample_rate, ipp):
                 return int(sample_rate*ipp*1e6 * r / 60) + int(sample_rate*ipp*1e6 * 1.2 / 60)
             def main(args):
                 experiment = args.experiment
                 fp = open(os.path.join(PATH, experiment, 'experiment.conf'))
                 conf = json.loads(fp.read())
                 ipp_km = conf['usrp_tx']['ipp']
                 ipp = ipp_km * 2 /300000
                 sample_rate  = conf['usrp_rx']['sample_rate']
                 axis = ['0' if x=='elevation' else '1' for x in conf['pedestal']['axis']]      # AZIMUTH 1 ELEVACION 0
                 speed_axis = conf['pedestal']['speed']
                 steps = conf['pedestal']['table']
                 time_offset = args.time_offset
                 parameters = args.parameters
                 start_date = experiment.split('@')[1].split('T')[0].replace('-', '/')
                 end_date = start_date
                 if args.start_time:
                     start_time = args.start_time
                 else:
                     start_time = experiment.split('@')[1].split('T')[1].replace('-', ':')
                 end_time = '23:59:59'
                 N = int(1/(speed_axis[0]*ipp))                                               # 1 GRADO DE RESOLUCION
                 path = os.path.join(PATH, experiment, 'rawdata')
                 path_ped = os.path.join(PATH, experiment, 'position')
                 path_plots = os.path.join(PATH, experiment, 'plotsC0N'+str(args.range))
                 path_save = os.path.join(PATH, experiment, 'paramC0N'+str(args.range))
                 RMIX = 1.62
                 H0 = -1.68
                 MASK = 0.3
                 from schainpy.controller import Project
                 project = Project()
                 project.setup(id='1', name='Sophy', description='sophy proc')
                 reader = project.addReadUnit(datatype='DigitalRFReader',
                     path=path,
                     startDate=start_date,
                     endDate=end_date,
                     startTime=start_time,
                     endTime=end_time,
                     delay=30,
                     online=args.online,
                     walk=1,
                     ippKm = ipp_km,
                     getByBlock = 1,
                     nProfileBlocks = N,
                 )
                 if not conf['usrp_tx']['enable_2']: # One Pulse
                     voltage = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
                     if conf['usrp_tx']['code_type_1'] != 'None':
                         codes = [ c.strip() for c in conf['usrp_tx']['code_1'].split(',')]
                         code = []
                         for c in codes:
                             code.append([int(x) for x in c])
                         op = voltage.addOperation(name='Decoder', optype='other')
                         op.addParameter(name='code', value=code)
                         op.addParameter(name='nCode', value=len(code), format='int')
                         op.addParameter(name='nBaud', value=len(code[0]), format='int')
                         op = voltage.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
                         op.addParameter(name='n', value=len(code), format='int')
                         ncode = len(code)
                     else:
                         ncode = 1
                         code = ['0']
                     op = voltage.addOperation(name='setH0')
                     op.addParameter(name='h0', value=H0)
                     if args.range > 0:
                         op = voltage.addOperation(name='selectHeights')
                         op.addParameter(name='minIndex', value='0', format='int')
                         op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
                     op = voltage.addOperation(name='PulsePair_vRF', optype='other')
                     op.addParameter(name='n', value=int(N)/ncode, format='int')
                     #op.addParameter(name='removeDC', value=1, format='int')
                     proc = project.addProcUnit(datatype='ParametersProc', inputId=voltage.getId())
                     opObj10 = proc.addOperation(name="WeatherRadar")
                     opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0]))
                     opObj10.addParameter(name='Pt',value=((1e-6/sample_rate)*len(code[0])/ipp)*200)
                     op = proc.addOperation(name='PedestalInformation')
                     op.addParameter(name='path', value=path_ped, format='str')
                     op.addParameter(name='interval', value='0.04')
                     op.addParameter(name='time_offset', value=time_offset)
                     op.addParameter(name='mode', value='PPI')
                     for param in parameters:
                         op = proc.addOperation(name='Block360')
                         op.addParameter(name='runNextOp', value=True)
                         op= proc.addOperation(name='WeatherParamsPlot')
                         if args.save: op.addParameter(name='save', value=path_plots, format='str')
                         op.addParameter(name='save_period', value=-1)
                         op.addParameter(name='show', value=args.show)
                         op.addParameter(name='channels', value='1,')
                         op.addParameter(name='zmin', value=PARAM[param]['zmin'])
                         op.addParameter(name='zmax', value=PARAM[param]['zmax'])
                         op.addParameter(name='attr_data', value=param, format='str')
                         op.addParameter(name='labels', value=[PARAM[param]['label']])
                         op.addParameter(name='save_code', value=param)
                         op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
                         op.addParameter(name='colormap', value=PARAM[param]['colormap'])
                         op.addParameter(name='bgcolor', value='black')
                         if MASK: op.addParameter(name='mask', value=MASK, format='float')
                         if args.server:
                             op.addParameter(name='server', value='0.0.0.0:4444')
                             op.addParameter(name='exp_code', value='400')
                         desc = {
                                 'Data': {
                                     param: PARAM[param]['wrname'],
                                     'utctime': 'time'
                                 },
                                  'Metadata': {
                                     'heightList': 'range',
                                     'data_azi': 'azimuth',
                                     'data_ele': 'elevation',
                                     'mode_op': 'scan_type',
                                     'h0': 'range_correction',
                                 }
                             }
                         if args.save:
                             opObj10 = proc.addOperation(name='HDFWriter')
                             writer.addParameter(name='path', value=path_save, format='str')
                             writer.addParameter(name='Reset', value=True)
                             writer.addParameter(name='setType', value='weather')
                             writer.addParameter(name='description', value=json.dumps(desc))
                             writer.addParameter(name='blocksPerFile', value='1',format='int')
                             writer.addParameter(name='metadataList', value='heightList,data_azi,data_ele,mode_op,latitude,longitude,altitude,heading,radar_name,institution,contact,h0,range_unit')
                             writer.addParameter(name='dataList', value='{},utctime'.format(param))
                             writer.addParameter(name='mask', value=MASK, format='float')
                             # meta
                             writer.addParameter(name='latitude', value='-12.040436')
                             writer.addParameter(name='longitude', value='-75.295893')
                             writer.addParameter(name='altitude', value='3379.2147')
                             writer.addParameter(name='heading', value='0')
                             writer.addParameter(name='radar_name', value='SOPHy')
                             writer.addParameter(name='institution', value='IGP')
                             writer.addParameter(name='contact', value='dscipion@igp.gob.pe')
                             writer.addParameter(name='created_by', value='Signal Chain (https://pypi.org/project/schainpy/)')
                             writer.addParameter(name='range_unit', value='km')
                 else: #Two pulses
                     voltage1 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
                     op = voltage1.addOperation(name='ProfileSelector')
                     op.addParameter(name='profileRangeList', value='0,{}'.format(conf['usrp_tx']['repetitions_1']-1))
                     if conf['usrp_tx']['code_type_1'] != 'None':
                         codes = [ c.strip() for c in conf['usrp_tx']['code_1'].split(',')]
                         code = []
                         for c in codes:
                             code.append([int(x) for x in c])
                         op = voltage1.addOperation(name='Decoder', optype='other')
                         op.addParameter(name='code', value=code)
                         op.addParameter(name='nCode', value=len(code), format='int')
                         op.addParameter(name='nBaud', value=len(code[0]), format='int')
                     else:
                         code = ['0']
                     op = voltage1.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
                     op.addParameter(name='n', value=2, format='int')
                     if args.range > 0:
                         op = voltage1.addOperation(name='selectHeights')
                         op.addParameter(name='minIndex', value='0', format='int')
                         op.addParameter(name='maxIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
                     op = voltage1.addOperation(name='setH0')
                     op.addParameter(name='h0', value=H0, format='float')
                     op = voltage1.addOperation(name='PulsePair_vRF', optype='other')
                     op.addParameter(name='n', value=int(conf['usrp_tx']['repetitions_1'])/2, format='int')
                     #op.addParameter(name='removeDC', value=1, format='int')
                     proc1 = project.addProcUnit(datatype='ParametersProc', inputId=voltage1.getId())
                     proc1.addParameter(name='runNextUnit', value=True)
                     opObj10 = proc1.addOperation(name="WeatherRadar")
                     opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0]))
-                    opObj10.addParameter(name='Pt',value=((1e-6/sample_rate)*len(code[0])/ipp)*200)
+                    opObj10.addParameter(name='Pt',value=200)
                     op = proc1.addOperation(name='PedestalInformation')
                     op.addParameter(name='path', value=path_ped, format='str')
                     op.addParameter(name='interval', value='0.04')
                     op.addParameter(name='time_offset', value=time_offset)
                     op.addParameter(name='mode', value='PPI')
                     op = proc1.addOperation(name='Block360')
                     op.addParameter(name='attr_data', value='data_param')
                     op.addParameter(name='runNextOp', value=True)
                     voltage2 = project.addProcUnit(datatype='VoltageProc', inputId=reader.getId())
                     op = voltage2.addOperation(name='ProfileSelector')
                     op.addParameter(name='profileRangeList', value='{},{}'.format(conf['usrp_tx']['repetitions_1'], conf['usrp_tx']['repetitions_1']+conf['usrp_tx']['repetitions_2']-1))
                     if conf['usrp_tx']['code_type_2']:
                         codes = [ c.strip() for c in conf['usrp_tx']['code_2'].split(',')]
                         code = []
                         for c in codes:
                             code.append([int(x) for x in c])
                         op = voltage2.addOperation(name='Decoder', optype='other')
                         op.addParameter(name='code', value=code)
                         op.addParameter(name='nCode', value=len(code), format='int')
                         op.addParameter(name='nBaud', value=len(code[0]), format='int')
                         op = voltage2.addOperation(name='CohInt', optype='other') #Minimo integrar 2 perfiles por ser codigo complementario
                         op.addParameter(name='n', value=len(code), format='int')
                         ncode = len(code)
                     else:
                         ncode = 1
                     if args.range > 0:
                         op = voltage2.addOperation(name='selectHeights')
                         op.addParameter(name='minIndex', value=max_index(RMIX, sample_rate, ipp), format='int')
                         op.addParameter(name='maxIndex', value=max_index(args.range, sample_rate, ipp), format='int')
                     op = voltage2.addOperation(name='setH0')
                     op.addParameter(name='h0', value=H0, format='float')
                     op = voltage2.addOperation(name='PulsePair_vRF', optype='other')
                     op.addParameter(name='n', value=int(conf['usrp_tx']['repetitions_2'])/ncode, format='int')
                     #op.addParameter(name='removeDC', value=1, format='int')
                     proc2 = project.addProcUnit(datatype='ParametersProc', inputId=voltage2.getId())
                     proc2.addParameter(name='runNextUnit', value=True)
                     opObj10 = proc2.addOperation(name="WeatherRadar")
                     opObj10.addParameter(name='tauW',value=(1e-6/sample_rate)*len(code[0]))
-                    opObj10.addParameter(name='Pt',value=((1e-6/sample_rate)*len(code[0])/ipp)*200)
+                    opObj10.addParameter(name='Pt',value=200)
                     op = proc2.addOperation(name='PedestalInformation')
                     op.addParameter(name='path', value=path_ped, format='str')
                     op.addParameter(name='interval', value='0.04')
                     op.addParameter(name='time_offset', value=time_offset)
                     op.addParameter(name='mode', value='PPI')
                     op = proc2.addOperation(name='Block360')
                     op.addParameter(name='attr_data', value='data_param')
                     op.addParameter(name='runNextOp', value=True)
                     merge = project.addProcUnit(datatype='MergeProc', inputId=[proc1.getId(), proc2.getId()])
                     merge.addParameter(name='attr_data', value='data_param')
                     merge.addParameter(name='mode', value='7') #RM
                     for param in parameters:
                         if args.plot:
                             op= merge.addOperation(name='WeatherParamsPlot')
                             if args.save:
                                 op.addParameter(name='save', value=path_plots, format='str')
                             op.addParameter(name='save_period', value=-1)
                             op.addParameter(name='show', value=args.show)
                             op.addParameter(name='channels', value='0,')
                             op.addParameter(name='zmin', value=PARAM[param]['zmin'], format='int')
                             op.addParameter(name='zmax', value=PARAM[param]['zmax'], format='int')
                             op.addParameter(name='attr_data', value=param, format='str')
                             op.addParameter(name='labels', value=[PARAM[param]['label']])
                             op.addParameter(name='save_code', value=param)
                             op.addParameter(name='cb_label', value=PARAM[param]['cb_label'])
                             op.addParameter(name='colormap', value=PARAM[param]['colormap'])
                             op.addParameter(name='bgcolor', value='black')
                             if MASK: op.addParameter(name='mask', value=MASK, format='float')
                             if args.server:
                                 op.addParameter(name='server', value='0.0.0.0:4444')
                                 op.addParameter(name='exp_code', value='400')
                         desc = {
                                 'Data': {
                                     'data_param': {PARAM[param]['wrname']: ['H', 'V']},
                                     'utctime': 'time'
                                 },
                                  'Metadata': {
                                     'heightList': 'range',
                                     'data_azi': 'azimuth',
                                     'data_ele': 'elevation',
                                     'mode_op': 'scan_type',
                                     'h0': 'range_correction',
                                 }
                             }
                         if args.save:
                             writer = merge.addOperation(name='HDFWriter')
                             writer.addParameter(name='path', value=path_save, format='str')
                             writer.addParameter(name='Reset', value=True)
                             writer.addParameter(name='setType', value='weather')
                             writer.addParameter(name='description', value=json.dumps(desc))
                             writer.addParameter(name='blocksPerFile', value='1',format='int')
                             writer.addParameter(name='metadataList', value='heightList,data_azi,data_ele,mode_op,latitude,longitude,altitude,heading,radar_name,institution,contact,h0,range_unit')
                             writer.addParameter(name='dataList', value='data_param,utctime')
                             writer.addParameter(name='weather_var', value=param)
                             writer.addParameter(name='mask', value=MASK, format='float')
                             # meta
                             writer.addParameter(name='latitude', value='-12.040436')
                             writer.addParameter(name='longitude', value='-75.295893')
                             writer.addParameter(name='altitude', value='3379.2147')
                             writer.addParameter(name='heading', value='0')
                             writer.addParameter(name='radar_name', value='SOPHy')
                             writer.addParameter(name='institution', value='IGP')
                             writer.addParameter(name='contact', value='dscipion@igp.gob.pe')
                             writer.addParameter(name='created_by', value='Signal Chain (https://pypi.org/project/schainpy/)')
                             writer.addParameter(name='range_unit', value='km')
                 project.start()
             if __name__ == '__main__':
                 parser = argparse.ArgumentParser(description='Script to process SOPHy data.')
                 parser.add_argument('experiment',
                                     help='Experiment name')
                 parser.add_argument('--parameters', nargs='*', default=['S'],
                                     help='Variables to process: P, Z, V')
                 parser.add_argument('--time_offset', default=0,
                                     help='Fix time offset')
                 parser.add_argument('--range', default=0, type=float,
                                     help='Max range to plot')
                 parser.add_argument('--save', action='store_true',
                                     help='Create output files')
                 parser.add_argument('--plot', action='store_true',
                                     help='Create plot files')
                 parser.add_argument('--show', action='store_true',
                                     help='Show matplotlib plot.')
                 parser.add_argument('--online', action='store_true',
                                     help='Set online mode.')
                 parser.add_argument('--server', action='store_true',
                                     help='Send to realtime')
                 parser.add_argument('--start_time', default='',
                                     help='Set start time.')
                 args = parser.parse_args()
                 main(args)
             # python sophy_proc.py  HYO_PM@2022-06-09T15-05-12  --parameters V --plot --save --show --range 36S

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