diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 62e57ca..df1f9b8 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -235,7 +235,6 @@ class JROData(GenericData): return delta def getltctime(self): - if self.useLocalTime: return self.utctime - self.timeZone * 60 @@ -1032,7 +1031,7 @@ class Correlation(JROData): normFactor = property(getNormFactor, "I'm the 'normFactor property'") -class Parameters(Spectra): +class Parameters(JROData): experimentInfo = None # Information about the experiment # Information from previous data @@ -1087,7 +1086,6 @@ class Parameters(Spectra): return datatime def getTimeInterval(self): - if hasattr(self, 'timeInterval1'): return self.timeInterval1 else: @@ -1210,7 +1208,7 @@ class PlotterData(object): ''' Update data object with new dataOut ''' - + print("UPDATE ESTOYU AQUI") self.profileIndex = dataOut.profileIndex self.tm = tm self.type = dataOut.type @@ -1278,6 +1276,14 @@ class PlotterData(object): buffer = dataOut.dataPP_WIDTH self.flagDataAsBlock = dataOut.flagDataAsBlock self.nProfiles = dataOut.nProfiles + if plot == 'pp_rti_power': + buffer = 10*numpy.log10(dataOut.dataPP_POWER) + if plot == 'pp_rti_signal': + buffer = 10*numpy.log10(dataOut.dataPP_POW) + if plot == 'weather': + print("ESTOY AQUI") + print("jrodata",dataOut.data_360.shape) + buffer = 10*numpy.log10(dataOut.data_360/(625**2)) if plot == 'spc': self.data['spc'][tm] = buffer diff --git a/schainpy/model/graphics/jroplot_base.py b/schainpy/model/graphics/jroplot_base.py index ccfbbc3..1696b11 100644 --- a/schainpy/model/graphics/jroplot_base.py +++ b/schainpy/model/graphics/jroplot_base.py @@ -9,6 +9,7 @@ from queue import Queue from functools import wraps from threading import Thread import matplotlib +import wradlib as wrl if 'BACKEND' in os.environ: matplotlib.use(os.environ['BACKEND']) @@ -233,7 +234,7 @@ class Plot(Operation): self.__throttle_plot = apply_throttle(self.throttle) self.data = PlotterData( self.CODE, self.throttle, self.exp_code, self.buffering, snr=self.showSNR) - + if self.plot_server: if not self.plot_server.startswith('tcp://'): self.plot_server = 'tcp://{}'.format(self.plot_server) @@ -250,8 +251,7 @@ class Plot(Operation): ''' self.setup() - - self.time_label = 'LT' if self.localtime else 'UTC' + self.time_label = 'LT' if self.localtime else 'UTC' if self.width is None: self.width = 8 @@ -273,6 +273,7 @@ class Plot(Operation): facecolor='w') self.figures.append(fig) for n in range(self.nplots): + #cgax,ax,paax = wrl.vis.create_cg(data_w,r=r,az=azi, proj='cg') ax = fig.add_subplot(self.nrows, self.ncols, n + 1, polar=self.polar) ax.tick_params(labelsize=8) @@ -381,10 +382,10 @@ class Plot(Operation): xmax += time.timezone else: xmax = self.xmax - + ymin = self.ymin if self.ymin else numpy.nanmin(self.y) ymax = self.ymax if self.ymax else numpy.nanmax(self.y) - + for n, ax in enumerate(self.axes): if ax.firsttime: @@ -398,14 +399,14 @@ class Plot(Operation): ystep = ystep/5 ystep = ystep/(10**digD) - else: + else: ystep = ((ymax + (10**(dig)))//10**(dig))*(10**(dig)) ystep = ystep/5 - + if self.xaxis is not 'time': - + dig = int(numpy.log10(xmax)) - + if dig <= 0: digD = len(str(xmax)) - 2 xdec = xmax*(10**digD) @@ -414,8 +415,8 @@ class Plot(Operation): xstep = ((xdec + (10**(dig)))//10**(dig))*(10**(dig)) xstep = xstep*0.5 xstep = xstep/(10**digD) - - else: + + else: xstep = ((xmax + (10**(dig)))//10**(dig))*(10**(dig)) xstep = xstep/5 @@ -495,14 +496,14 @@ class Plot(Operation): self.plot() self.format() - + for n, fig in enumerate(self.figures): if self.nrows == 0 or self.nplots == 0: log.warning('No data', self.name) fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center') fig.canvas.manager.set_window_title(self.CODE) continue - + fig.canvas.manager.set_window_title('{} - {}'.format(self.title, self.getDateTime(self.data.max_time).strftime('%Y/%m/%d'))) fig.canvas.draw() @@ -512,7 +513,7 @@ class Plot(Operation): if self.save: self.save_figure(n) - + if self.plot_server: self.send_to_server() @@ -539,7 +540,7 @@ class Plot(Operation): figname = os.path.join( self.save, labels[n], - '{}_{}.png'.format( + '{}_{}.png'.format( labels[n], self.getDateTime(self.data.max_time).strftime( '%Y%m%d_%H%M%S' @@ -572,7 +573,7 @@ class Plot(Operation): return self.sender_time = self.data.tm - + attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax'] for attr in attrs: value = getattr(self, attr) @@ -589,7 +590,7 @@ class Plot(Operation): self.data.meta['interval'] = int(interval) # msg = self.data.jsonify(self.data.tm, self.plot_name, self.plot_type) self.sender_queue.put(self.data.tm) - + while True: if self.sender_queue.empty(): break @@ -628,7 +629,7 @@ class Plot(Operation): self.ncols: number of cols self.nplots: number of plots (channels or pairs) self.ylabel: label for Y axes - self.titles: list of axes title + self.titles: list of axes title ''' raise NotImplementedError @@ -638,7 +639,7 @@ class Plot(Operation): Must be defined in the child class ''' raise NotImplementedError - + def run(self, dataOut, **kwargs): ''' Main plotting routine @@ -646,7 +647,7 @@ class Plot(Operation): if self.isConfig is False: self.__setup(**kwargs) - + t = getattr(dataOut, self.attr_time) if dataOut.useLocalTime: @@ -657,7 +658,7 @@ class Plot(Operation): self.getDateTime = datetime.datetime.utcfromtimestamp if self.localtime: t -= time.timezone - + if self.xmin is None: self.tmin = t if 'buffer' in self.plot_type: @@ -665,8 +666,8 @@ class Plot(Operation): else: self.tmin = ( self.getDateTime(t).replace( - hour=int(self.xmin), - minute=0, + hour=int(self.xmin), + minute=0, second=0) - self.getDateTime(0)).total_seconds() self.data.setup() @@ -684,8 +685,8 @@ class Plot(Operation): tm -= time.timezone if dataOut.useLocalTime and not self.localtime: tm += time.timezone - - if self.data and (tm - self.tmin) >= self.xrange*60*60: + + if self.data and (tm - self.tmin) >= self.xrange*60*60: self.save_counter = self.save_period self.__plot() if 'time' in self.xaxis: @@ -697,7 +698,7 @@ class Plot(Operation): self.clear_figures() self.data.update(dataOut, tm) - + print("plotbase--",self.data['weather'].shape) if self.isPlotConfig is False: self.__setup_plot() self.isPlotConfig = True @@ -714,4 +715,3 @@ class Plot(Operation): self.__plot() if self.data and not self.data.flagNoData and self.pause: figpause(10) - diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index e275f12..7c1c660 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -6,6 +6,10 @@ from schainpy.model.graphics.jroplot_base import Plot, plt from schainpy.model.graphics.jroplot_spectra import SpectraPlot, RTIPlot, CoherencePlot from schainpy.utils import log +'''-------''' +import wradlib as wrl +'''-------''' + EARTH_RADIUS = 6.3710e3 @@ -69,6 +73,24 @@ class PowerPlot(RTIPlot): colormap = 'jet' plot_name = 'TotalPower' +class PPPowerPlot(RTIPlot): + ''' + Plot for Pulse Pair Power Data P =S +N + ''' + + CODE = 'pp_rti_power' + colormap = 'jet' + plot_name = 'TotalPP_Power' + +class PPSignalPlot(RTIPlot): + ''' + Plot for Pulse Pair Power Data S = P- N (0 moment) + ''' + + CODE = 'pp_rti_signal' + colormap = 'jet' + plot_name = 'TotalPP_Signal' + class SpectralWidthPlot(RTIPlot): ''' @@ -339,3 +361,53 @@ class PolarMapPlot(Plot): self.titles = ['{} {}'.format( self.data.parameters[x], title) for x in self.channels] + +class WeatherPlot(Plot): + CODE = 'weather' + plot_name = 'Weather' + plot_type = 'ppistyle' + + def setup(self): + #self.xaxis = 'Range (Km)' + self.ncols = 1 + self.nrows = 1 + self.nplots = 1 + self.ylabel = 'Range (Km)' + self.titles = ['Weather'] + self.colorbar = False + self.width = 8 + self.height = 8 + self.ini = 50 + + def plot(self): + stoprange = float(33*1.5) + rangestep = float(0.15) + r = numpy.arange(0, stoprange, rangestep) + ''' + self.x= r + self.y= r + ini = self.ini + self.ini = self.ini+100 + azi = numpy.linspace(0,359,360) + azi = azi + ini + azi = azi % 36 + ''' + azi = self.data['weather'][1] + data = self.data['weather'][0] + data_w= data[:,-1,:] + try: + zeros=numpy.zeros([(360-data_w.shape[0]),data_w.shape[1]]) + print("zeros",zeros.shape) + data_w= numpy.vstack([data_w,zeros]) + except: + pass + + for i,ax in enumerate(self.axes): + if ax.firsttime: + plt.clf() + cgax, pm = wrl.vis.plot_ppi(data_w,r=r,az=azi,fig=self.figures[0], proj='cg') + else: + plt.clf() + cgax, pm = wrl.vis.plot_ppi(data_w,r=r,az=azi,fig=self.figures[0], proj='cg') + + print("cgax",cgax) diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 3ef5c6e..724b56d 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -121,7 +121,7 @@ class CrossSpectraPlot(Plot): else: x = self.data.xrange[2] self.xlabel = "Velocity (m/s)" - + self.titles = [] y = self.data.heights @@ -134,17 +134,17 @@ class CrossSpectraPlot(Plot): noise = self.data['noise'][:,-1] pair = self.data.pairs[n] ax = self.axes[4 * n] - 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(nspc) - self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc) + 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(nspc) + self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc) ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T, vmin=self.zmin, vmax=self.zmax, cmap=plt.get_cmap(self.colormap) - ) - else: + ) + else: ax.plt.set_array(nspc[pair[0]].T.ravel()) self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]])) @@ -155,10 +155,10 @@ class CrossSpectraPlot(Plot): vmax=self.zmax, cmap=plt.get_cmap(self.colormap) ) - else: + else: ax.plt.set_array(nspc[pair[1]].T.ravel()) self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]])) - + 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 @@ -174,13 +174,13 @@ class CrossSpectraPlot(Plot): ax.plt.set_array(coh.T.ravel()) self.titles.append( 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1])) - + ax = self.axes[4 * n + 3] if ax.firsttime: ax.plt = ax.pcolormesh(x, y, phase.T, vmin=-180, vmax=180, - cmap=plt.get_cmap(self.colormap_phase) + cmap=plt.get_cmap(self.colormap_phase) ) else: ax.plt.set_array(phase.T.ravel()) @@ -285,7 +285,7 @@ class PhasePlot(CoherencePlot): class NoisePlot(Plot): ''' - Plot for noise + Plot for noise ''' CODE = 'noise' @@ -349,10 +349,10 @@ class PowerProfilePlot(Plot): self.y = y x = self.data['spcprofile'] - + 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)) @@ -506,7 +506,7 @@ class BeaconPhase(Plot): server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): - if dataOut.flagNoData: + if dataOut.flagNoData: return dataOut if not isTimeInHourRange(dataOut.datatime, xmin, xmax): @@ -648,4 +648,4 @@ class BeaconPhase(Plot): thisDatetime=thisDatetime, update_figfile=update_figfile) - return dataOut \ No newline at end of file + return dataOut diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index 0e176ba..0d9a457 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -11,6 +11,8 @@ from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecora from schainpy.model.io.jroIO_base import * from schainpy.utils import log +import wradlib as wrl + class ParamReader(JRODataReader,ProcessingUnit): ''' @@ -183,17 +185,18 @@ class ParamReader(JRODataReader,ProcessingUnit): except IOError: traceback.print_exc() raise IOError("The file %s can't be opened" %(filename)) - + #In case has utctime attribute grp2 = grp1['utctime'] # thisUtcTime = grp2.value[0] - 5*3600 #To convert to local time thisUtcTime = grp2.value[0] - + #thisUtcTime = grp2 fp.close() if self.timezone == 'lt': thisUtcTime -= 5*3600 + thisDatetime = datetime.datetime.fromtimestamp(thisUtcTime[0] + 5*3600) thisDate = thisDatetime.date() thisTime = thisDatetime.time() @@ -266,8 +269,9 @@ class ParamReader(JRODataReader,ProcessingUnit): endTime = self.endTime grp = fp['Data'] + #grp2 = grp['utctime'][()] # nuevo + #thisUtcTime = grp2 # nuevo thisUtcTime = grp['utctime'].value.astype(numpy.float)[0] - #ERROOOOR if self.timezone == 'lt': thisUtcTime -= 5*3600 @@ -304,14 +308,14 @@ class ParamReader(JRODataReader,ProcessingUnit): listMetaname = [] listMetadata = [] + listShapes = {} for item in list(gp.items()): name = item[0] - if name=='array dimensions': table = gp[name][:] listShapes = {} for shapes in table: - listShapes[shapes[0]] = numpy.array([shapes[1],shapes[2],shapes[3],shapes[4],shapes[5]]) + listShapes[shapes[0].decode()] = numpy.array([shapes[1],shapes[2],shapes[3],shapes[4],shapes[5]]) else: data = gp[name].value listMetaname.append(name) @@ -328,11 +332,9 @@ class ParamReader(JRODataReader,ProcessingUnit): grp = self.fp['Data'] listdataname = [] listdata = [] - for item in list(grp.items()): name = item[0] listdataname.append(name) - array = self.__setDataArray(grp[name],self.listShapes[name]) listdata.append(array) @@ -341,7 +343,6 @@ class ParamReader(JRODataReader,ProcessingUnit): return def __setDataArray(self, dataset, shapes): - nDims = shapes[0] nDim2 = shapes[1] #Dimension 0 nDim1 = shapes[2] #Dimension 1, number of Points or Parameters @@ -407,7 +408,6 @@ class ParamReader(JRODataReader,ProcessingUnit): listShapes = self.listShapes blockIndex = self.blockIndex - # blockList = self.blockList for i in range(len(listMeta)): setattr(self.dataOut,listMetaname[i],listMeta[i]) @@ -497,7 +497,7 @@ class ParamWriter(Operation): setType = None def __init__(self): - + Operation.__init__(self) return @@ -530,9 +530,9 @@ class ParamWriter(Operation): dsDict['variable'] = self.dataList[i] #--------------------- Conditionals ------------------------ #There is no data - + if dataAux is None: - + return 0 if isinstance(dataAux, (int, float, numpy.integer, numpy.float)): @@ -704,7 +704,7 @@ class ParamWriter(Operation): return False def setNextFile(self): - + ext = self.ext path = self.path setFile = self.setFile @@ -785,7 +785,7 @@ class ParamWriter(Operation): for j in range(dsInfo['dsNumber']): dsInfo = dsList[i] tableName = dsInfo['dsName'] - + if dsInfo['nDim'] == 3: shape = dsInfo['shape'].astype(int) @@ -954,7 +954,7 @@ class ParamWriter(Operation): self.dataOut = dataOut if not(self.isConfig): - self.setup(dataOut, path=path, blocksPerFile=blocksPerFile, + self.setup(dataOut, path=path, blocksPerFile=blocksPerFile, metadataList=metadataList, dataList=dataList, mode=mode, setType=setType) @@ -963,7 +963,7 @@ class ParamWriter(Operation): self.putData() return - + class ParameterReader(Reader, ProcessingUnit): @@ -992,55 +992,61 @@ class ParameterReader(Reader, ProcessingUnit): self.set_kwargs(**kwargs) if not self.ext.startswith('.'): - self.ext = '.{}'.format(self.ext) + self.ext = '.{}'.format(self.ext) if self.online: log.log("Searching files in online mode...", self.name) for nTries in range(self.nTries): fullpath = self.searchFilesOnLine(self.path, self.startDate, - self.endDate, self.expLabel, self.ext, self.walk, + self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt) try: fullpath = next(fullpath) except: fullpath = None - + if fullpath: break log.warning( 'Waiting {} sec for a valid file in {}: try {} ...'.format( - self.delay, self.path, nTries + 1), + self.delay, self.path, nTries + 1), self.name) time.sleep(self.delay) if not(fullpath): raise schainpy.admin.SchainError( - 'There isn\'t any valid file in {}'.format(self.path)) + 'There isn\'t any valid file in {}'.format(self.path)) pathname, filename = os.path.split(fullpath) self.year = int(filename[1:5]) self.doy = int(filename[5:8]) - self.set = int(filename[8:11]) - 1 + self.set = int(filename[8:11]) - 1 else: log.log("Searching files in {}".format(self.path), self.name) - self.filenameList = self.searchFilesOffLine(self.path, self.startDate, + self.filenameList = self.searchFilesOffLine(self.path, self.startDate, self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt) - - self.setNextFile() + self.setNextFile() + ''' + # try to implement wradlib + pathname,filename = os.path.split(self.filename) + f = wrl.util.get_wradlib_data_file(self.filename) + fcontent = wrl.io.read_generic_hdf5(f) + print(pathname,filename) + ''' return def readFirstHeader(self): '''Read metadata and data''' - self.__readMetadata() + self.__readMetadata() self.__readData() self.__setBlockList() self.blockIndex = 0 - + return def __setBlockList(self): @@ -1057,19 +1063,27 @@ class ParameterReader(Reader, ProcessingUnit): startTime = self.startTime endTime = self.endTime - index = self.listDataname.index('utctime') thisUtcTime = self.listData[index] self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1]) - + try: + index = self.listMetaname.index('timeZone') + tz = self.listMeta[index] + if tz == 0: + print("[Metadata] Attribute timeZone = %d, dataOut in utc"%(tz)) + else: + print("[Metadata] Attribute timeZone = %d, dataOut in LT"%(tz)) + thisUtcTime += tz*60 + except: + print("[Metadata] There is no attribute timeZone , by default assume dataOut was recorded in LT") + tz=300 + thisUtcTime += tz*60 if self.timezone == 'lt': thisUtcTime -= 5*3600 - thisDatetime = datetime.datetime.fromtimestamp(thisUtcTime[0] + 5*3600) - + thisDatetime = datetime.datetime.fromtimestamp(thisUtcTime[0]) thisDate = thisDatetime.date() thisTime = thisDatetime.time() - startUtcTime = (datetime.datetime.combine(thisDate,startTime) - datetime.datetime(1970, 1, 1)).total_seconds() endUtcTime = (datetime.datetime.combine(thisDate,endTime) - datetime.datetime(1970, 1, 1)).total_seconds() @@ -1099,7 +1113,7 @@ class ParameterReader(Reader, ProcessingUnit): else: data = gp[name].value listMetaname.append(name) - listMetadata.append(data) + listMetadata.append(data) elif self.metadata: metadata = json.loads(self.metadata) listShapes = {} @@ -1115,7 +1129,7 @@ class ParameterReader(Reader, ProcessingUnit): self.listShapes = listShapes self.listMetaname = listMetaname - self.listMeta = listMetadata + self.listMeta = listMetadata return @@ -1123,7 +1137,7 @@ class ParameterReader(Reader, ProcessingUnit): listdataname = [] listdata = [] - + if 'Data' in self.fp: grp = self.fp['Data'] for item in list(grp.items()): @@ -1137,7 +1151,7 @@ class ParameterReader(Reader, ProcessingUnit): for i in range(dim): array.append(grp[name]['table{:02d}'.format(i)].value) array = numpy.array(array) - + listdata.append(array) elif self.metadata: metadata = json.loads(self.metadata) @@ -1160,9 +1174,8 @@ class ParameterReader(Reader, ProcessingUnit): self.listDataname = listdataname self.listData = listdata return - - def getData(self): + def getData(self): for i in range(len(self.listMeta)): setattr(self.dataOut, self.listMetaname[i], self.listMeta[i]) @@ -1174,9 +1187,16 @@ class ParameterReader(Reader, ProcessingUnit): setattr(self.dataOut, self.listDataname[j], self.listData[j][:,self.blockIndex]) self.dataOut.paramInterval = self.interval + + try: + thisDatetime = datetime.datetime.fromtimestamp(self.dataOut.utctime) + '''REMEMBER THIS''' + #print("[Reading] Block No. %d -> %s" % (self.blockIndex,thisDatetime.ctime())) + except: + pass + self.dataOut.flagNoData = False self.blockIndex += 1 - return def run(self, **kwargs): @@ -1230,7 +1250,7 @@ class ParameterWriter(Operation): lastTime = None def __init__(self): - + Operation.__init__(self) return @@ -1257,7 +1277,7 @@ class ParameterWriter(Operation): dsDict['nDim'] = len(dataAux.shape) dsDict['shape'] = dataAux.shape dsDict['dsNumber'] = dataAux.shape[0] - + dsList.append(dsDict) tableList.append((self.dataList[i], dsDict['nDim'])) @@ -1274,7 +1294,7 @@ class ParameterWriter(Operation): self.lastTime = currentTime self.currentDay = dataDay return False - + timeDiff = currentTime - self.lastTime #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora @@ -1292,7 +1312,7 @@ class ParameterWriter(Operation): self.dataOut = dataOut if not(self.isConfig): - self.setup(path=path, blocksPerFile=blocksPerFile, + self.setup(path=path, blocksPerFile=blocksPerFile, metadataList=metadataList, dataList=dataList, setType=setType) @@ -1301,9 +1321,9 @@ class ParameterWriter(Operation): self.putData() return - + def setNextFile(self): - + ext = self.ext path = self.path setFile = self.setFile @@ -1369,17 +1389,17 @@ class ParameterWriter(Operation): return def writeData(self, fp): - + grp = fp.create_group("Data") dtsets = [] data = [] - + for dsInfo in self.dsList: if dsInfo['nDim'] == 0: ds = grp.create_dataset( - dsInfo['variable'], + dsInfo['variable'], (self.blocksPerFile, ), - chunks=True, + chunks=True, dtype=numpy.float64) dtsets.append(ds) data.append((dsInfo['variable'], -1)) @@ -1387,7 +1407,7 @@ class ParameterWriter(Operation): sgrp = grp.create_group(dsInfo['variable']) for i in range(dsInfo['dsNumber']): ds = sgrp.create_dataset( - 'table{:02d}'.format(i), + 'table{:02d}'.format(i), (self.blocksPerFile, ) + dsInfo['shape'][1:], chunks=True) dtsets.append(ds) @@ -1395,7 +1415,7 @@ class ParameterWriter(Operation): fp.flush() log.log('Creating file: {}'.format(fp.filename), self.name) - + self.ds = dtsets self.data = data self.firsttime = True diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index f521bb6..bc9f4d9 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -109,6 +109,13 @@ class ParametersProc(ProcessingUnit): self.dataOut.flagNoData = False self.dataOut.utctimeInit = self.dataIn.utctime self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds + + if hasattr(self.dataIn, 'flagDataAsBlock'): + self.dataOut.flagDataAsBlock = self.dataIn.flagDataAsBlock + + if hasattr(self.dataIn, 'profileIndex'): + self.dataOut.profileIndex = self.dataIn.profileIndex + if hasattr(self.dataIn, 'dataPP_POW'): self.dataOut.dataPP_POW = self.dataIn.dataPP_POW @@ -183,7 +190,6 @@ class ParametersProc(ProcessingUnit): if self.dataIn.type == "Parameters": self.dataOut.copy(self.dataIn) self.dataOut.flagNoData = False - return True self.__updateObjFromInput() @@ -592,7 +598,86 @@ class GaussianFit(Operation): def misfit2(self,state,y_data,x,num_intg): return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.) +class WeatherRadar(Operation): + ''' + Function tat implements Weather Radar operations- + Input: + Output: + Parameters affected: + ''' + isConfig = False + + def __init__(self): + Operation.__init__(self) + def setup(self,dataOut,Pt=0,Gt=0,Gr=0,lambda_=0, aL=0, + tauW= 0,thetaT=0,thetaR=0,Km =0): + self.nCh = dataOut.nChannels + self.nHeis = dataOut.nHeights + deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] + self.Range = numpy.arange(dataOut.nHeights)*deltaHeight + dataOut.heightList[0] + self.Range = self.Range.reshape(1,self.nHeis) + self.Range = numpy.tile(self.Range,[self.nCh,1]) + '''-----------1 Constante del Radar----------''' + self.Pt = Pt + self.Gt = Gt + self.Gr = Gr + self.lambda_ = lambda_ + self.aL = aL + self.tauW = tauW + self.thetaT = thetaT + self.thetaR = thetaR + self.Km = Km + Numerator = ((4*numpy.pi)**3 * aL**2 * 16 *numpy.log(2)) + Denominator = (Pt * Gt * Gr * lambda_**2 * SPEED_OF_LIGHT * tauW * numpy.pi*thetaT*thetaR) + self.RadarConstant = Numerator/Denominator + '''-----------2 Reflectividad del Radar y Factor de Reflectividad------''' + self.n_radar = numpy.zeros((self.nCh,self.nHeis)) + self.Z_radar = numpy.zeros((self.nCh,self.nHeis)) + + def setMoments(self,dataOut,i): + + type = dataOut.inputUnit + nCh = dataOut.nChannels + nHeis= dataOut.nHeights + data_param = numpy.zeros((nCh,4,nHeis)) + if type == "Voltage": + data_param[:,0,:] = dataOut.dataPP_POW/(dataOut.nCohInt**2) + data_param[:,1,:] = dataOut.dataPP_DOP + data_param[:,2,:] = dataOut.dataPP_WIDTH + data_param[:,3,:] = dataOut.dataPP_SNR + if type == "Spectra": + data_param[:,0,:] = dataOut.data_POW + data_param[:,1,:] = dataOut.data_DOP + data_param[:,2,:] = dataOut.data_WIDTH + data_param[:,3,:] = dataOut.data_SNR + + return data_param[:,i,:] + + + def run(self,dataOut,Pt=25,Gt=200.0,Gr=50.0,lambda_=0.32, aL=2.5118, + tauW= 4.0e-6,thetaT=0.165,thetaR=0.367,Km =0.93): + + if not self.isConfig: + self.setup(dataOut= dataOut,Pt=25,Gt=200.0,Gr=50.0,lambda_=0.32, aL=2.5118, + tauW= 4.0e-6,thetaT=0.165,thetaR=0.367,Km =0.93) + self.isConfig = True + '''-----------------------------Potencia de Radar -Signal S-----------------------------''' + Pr = self.setMoments(dataOut,0) + + for R in range(self.nHeis): + self.n_radar[:,R] = self.RadarConstant*Pr[:,R]* (self.Range[:,R])**2 + + self.Z_radar[:,R] = self.n_radar[:,R]* self.lambda_**4/( numpy.pi**5 * self.Km**2) + + '''----------- Factor de Reflectividad Equivalente lamda_ < 10 cm , lamda_= 3.2cm-------''' + Zeh = self.Z_radar + dBZeh = 10*numpy.log10(Zeh) + dataOut.factor_Zeh= dBZeh + self.n_radar = numpy.zeros((self.nCh,self.nHeis)) + self.Z_radar = numpy.zeros((self.nCh,self.nHeis)) + + return dataOut class PrecipitationProc(Operation): @@ -3998,3 +4083,249 @@ class SMOperations(): # error[indInvalid1] = 13 # # return heights, error + +class PulsePairVoltage(Operation): + ''' + Function PulsePair(Signal Power, Velocity) + The real component of Lag[0] provides Intensity Information + The imag component of Lag[1] Phase provides Velocity Information + + Configuration Parameters: + nPRF = Number of Several PRF + theta = Degree Azimuth angel Boundaries + + Input: + self.dataOut + lag[N] + Affected: + self.dataOut.spc + ''' + isConfig = False + __profIndex = 0 + __initime = None + __lastdatatime = None + __buffer = None + noise = None + __dataReady = False + n = None + __nch = 0 + __nHeis = 0 + removeDC = False + ipp = None + lambda_ = 0 + + def __init__(self,**kwargs): + Operation.__init__(self,**kwargs) + + def setup(self, dataOut, n = None, removeDC=False): + ''' + n= Numero de PRF's de entrada + ''' + self.__initime = None + self.__lastdatatime = 0 + self.__dataReady = False + self.__buffer = 0 + self.__profIndex = 0 + self.noise = None + self.__nch = dataOut.nChannels + self.__nHeis = dataOut.nHeights + self.removeDC = removeDC + self.lambda_ = 3.0e8/(9345.0e6) + self.ippSec = dataOut.ippSeconds + self.nCohInt = dataOut.nCohInt + print("IPPseconds",dataOut.ippSeconds) + + print("ELVALOR DE n es:", n) + if n == None: + raise ValueError("n should be specified.") + + if n != None: + if n<2: + raise ValueError("n should be greater than 2") + + self.n = n + self.__nProf = n + + self.__buffer = numpy.zeros((dataOut.nChannels, + n, + dataOut.nHeights), + dtype='complex') + + def putData(self,data): + ''' + Add a profile to he __buffer and increase in one the __profiel Index + ''' + self.__buffer[:,self.__profIndex,:]= data + self.__profIndex += 1 + return + + def pushData(self,dataOut): + ''' + Return the PULSEPAIR and the profiles used in the operation + Affected : self.__profileIndex + ''' + #················· Remove DC··································· + if self.removeDC==True: + mean = numpy.mean(self.__buffer,1) + tmp = mean.reshape(self.__nch,1,self.__nHeis) + dc= numpy.tile(tmp,[1,self.__nProf,1]) + self.__buffer = self.__buffer - dc + #··················Calculo de Potencia ························ + pair0 = self.__buffer*numpy.conj(self.__buffer) + pair0 = pair0.real + lag_0 = numpy.sum(pair0,1) + #··················Calculo de Ruido x canal···················· + self.noise = numpy.zeros(self.__nch) + for i in range(self.__nch): + daux = numpy.sort(pair0[i,:,:],axis= None) + self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt) + + self.noise = self.noise.reshape(self.__nch,1) + self.noise = numpy.tile(self.noise,[1,self.__nHeis]) + noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis) + noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1]) + #·················· Potencia recibida= P , Potencia senal = S , Ruido= N·· + #·················· P= S+N ,P=lag_0/N ································· + #···················· Power ·················································· + data_power = lag_0/(self.n*self.nCohInt) + #------------------ Senal ··················································· + data_intensity = pair0 - noise_buffer + data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt) + #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt) + for i in range(self.__nch): + for j in range(self.__nHeis): + if data_intensity[i][j] < 0: + data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j])) + + #················· Calculo de Frecuencia y Velocidad doppler········ + pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:]) + lag_1 = numpy.sum(pair1,1) + data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1) + data_velocity = (self.lambda_/2.0)*data_freq + + #················ Potencia promedio estimada de la Senal··········· + lag_0 = lag_0/self.n + S = lag_0-self.noise + + #················ Frecuencia Doppler promedio ····················· + lag_1 = lag_1/(self.n-1) + R1 = numpy.abs(lag_1) + + #················ Calculo del SNR·································· + data_snrPP = S/self.noise + for i in range(self.__nch): + for j in range(self.__nHeis): + if data_snrPP[i][j] < 1.e-20: + data_snrPP[i][j] = 1.e-20 + + #················· Calculo del ancho espectral ······················ + L = S/R1 + L = numpy.where(L<0,1,L) + L = numpy.log(L) + tmp = numpy.sqrt(numpy.absolute(L)) + data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L) + n = self.__profIndex +class Block360(Operation): + ''' + ''' + isConfig = False + __profIndex = 0 + __initime = None + __lastdatatime = None + __buffer = None + __dataReady = False + n = None + __nch = 0 + __nHeis = 0 + + def __init__(self,**kwargs): + Operation.__init__(self,**kwargs) + + def setup(self, dataOut, n = None): + ''' + n= Numero de PRF's de entrada + ''' + self.__initime = None + self.__lastdatatime = 0 + self.__dataReady = False + self.__buffer = 0 + self.__profIndex = 0 + + + print("ELVALOR DE n es:", n) + if n == None: + raise ValueError("n should be specified.") + + if n != None: + if n<2: + raise ValueError("n should be greater than 2") + + self.n = n + + self.__buffer = numpy.zeros(( n, dataOut.nHeights)) + + def putData(self,data): + ''' + Add a profile to he __buffer and increase in one the __profiel Index + ''' + self.__buffer[self.__profIndex,:]= data + self.__profIndex += 1 + return #················· Remove DC··································· + + def pushData(self,dataOut): + ''' + Return the PULSEPAIR and the profiles used in the operation + Affected : self.__profileIndex + ''' + + + data_360 = self.__buffer + n = self.__profIndex + + self.__buffer = numpy.zeros(( self.n,330)) + self.__profIndex = 0 + return data_360,n + + + def byProfiles(self,dataOut): + + self.__dataReady = False + data_360 = None + self.putData(data=dataOut.data_POW[0]) + if self.__profIndex == self.n: + data_360, n = self.pushData(dataOut=dataOut) + self.__dataReady = True + + return data_360 + + + def blockOp(self, dataOut, datatime= None): + + if self.__initime == None: + self.__initime = datatime + data_360 = self.byProfiles(dataOut) + self.__lastdatatime = datatime + + if data_360 is None: + return None, None + + avgdatatime = self.__initime + deltatime = datatime - self.__lastdatatime + self.__initime = datatime + + return data_360, avgdatatime + + def run(self, dataOut,n = None,**kwargs): + + if not self.isConfig: + self.setup(dataOut = dataOut, n = n , **kwargs) + self.isConfig = True + data_360, avgdatatime = self.blockOp(dataOut, dataOut.utctime) + dataOut.flagNoData = True + + if self.__dataReady: + dataOut.data_360 = data_360 # S + print("data_360",data_360.shape,avgdatatime) + dataOut.utctime = avgdatatime + dataOut.flagNoData = False + return dataOut diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index a54b500..4ba9971 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -1484,6 +1484,7 @@ class PulsePairVoltage(Operation): dataOut.dataPP_SNR = data_snrPP dataOut.dataPP_WIDTH = data_specwidth dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo. + dataOut.nProfiles = int(dataOut.nProfiles/n) dataOut.utctime = avgdatatime dataOut.flagNoData = False return dataOut diff --git a/schainpy/scripts/pedestal_client.py b/schainpy/scripts/pedestal_client.py index ec88d9e..50b376d 100644 --- a/schainpy/scripts/pedestal_client.py +++ b/schainpy/scripts/pedestal_client.py @@ -47,7 +47,7 @@ while(True): elev.append(ang_elev) time0.append(seconds) count +=1 - if count == 100: + if count == 360: timetuple=time.localtime() epoc = time.mktime(timetuple) #print(epoc) diff --git a/schainpy/scripts/test_sim00010.py b/schainpy/scripts/test_sim00010.py index 8710965..929f691 100644 --- a/schainpy/scripts/test_sim00010.py +++ b/schainpy/scripts/test_sim00010.py @@ -76,7 +76,7 @@ opObj10 = procUnitConfObjC.addOperation(name='ParameterWriter') opObj10.addParameter(name='path',value=figpath) #opObj10.addParameter(name='mode',value=0) opObj10.addParameter(name='blocksPerFile',value='100',format='int') -opObj10.addParameter(name='metadataList',value='utctimeInit,timeInterval',format='list') -opObj10.addParameter(name='dataList',value='data_POW,data_DOP,data_WIDTH,data_SNR')#,format='list' +opObj10.addParameter(name='metadataList',value='utctimeInit,paramInterval,heightList',format='list') +opObj10.addParameter(name='dataList',value='data_POW,data_DOP,utctime')#,format='list' controllerObj.start() diff --git a/schainpy/scripts/test_sim0009.py b/schainpy/scripts/test_sim0009.py index 43d6640..013b283 100644 --- a/schainpy/scripts/test_sim0009.py +++ b/schainpy/scripts/test_sim0009.py @@ -2,7 +2,7 @@ import os,sys import datetime import time from schainpy.controller import Project -path = '/home/alex/Downloads/NEW_WR2/spc16removeDC' +path = '/home/alex/Downloads/hdf5_testPP2' figpath = path desc = "Simulator Test" @@ -51,7 +51,7 @@ opObj11 = procUnitConfObjA.addOperation(name='PulsePairVoltage', optype='other') opObj11.addParameter(name='n', value='625', format='int')#10 opObj11.addParameter(name='removeDC', value=1, format='int') -#opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other') +#opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')#,dataPP_SNR,dataPP_WIDTH #opObj11 = procUnitConfObjA.addOperation(name='PulsepairSignalPlot', optype='other') @@ -61,13 +61,11 @@ opObj11.addParameter(name='removeDC', value=1, format='int') #opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other') procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId()) - - opObj10 = procUnitConfObjB.addOperation(name='ParameterWriter') opObj10.addParameter(name='path',value=figpath) -#opObj10.addParameter(name='mode',value=0) -opObj10.addParameter(name='blocksPerFile',value='100',format='int') -opObj10.addParameter(name='metadataList',value='utctimeInit,timeInterval',format='list') -opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,dataPP_SNR,dataPP_WIDTH')#,format='list' +#opObj10.addParameter(name='mode',value=2) +opObj10.addParameter(name='blocksPerFile',value='50',format='int') +opObj10.addParameter(name='metadataList',value='utctimeInit,paramInterval,heightList,profileIndex,flagDataAsBlock',format='list') +opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,utctime',format='list')#,format='list' controllerObj.start() diff --git a/schainpy/scripts/wr_integrador.py b/schainpy/scripts/wr_integrador.py index d719cb3..51edaa3 100644 --- a/schainpy/scripts/wr_integrador.py +++ b/schainpy/scripts/wr_integrador.py @@ -82,9 +82,9 @@ num_perfiles = int(var_tiempo/IPP) dir_pedestal = "/home/alex/Downloads/pedestal" #·········· DATA ADQ······························ if mode=="T": - dir_adq = "/home/alex/Downloads/hdf5_testPP/d2020194" # Time domain + dir_adq = "/home/alex/Downloads/hdf5_testPP/d2020204" # Time domain else: - dir_adq = "/home/alex/Downloads/hdf5_test/d2020194" # Frequency domain + dir_adq = "/home/alex/Downloads/hdf5_test/d2020204" # Frequency domain print( "Velocidad angular :", w) print( "Resolucion minima en grados :", alfa) @@ -103,7 +103,7 @@ print("utc_pedestal :",utc_pedestal) print("utc_adq :",utc_adq) #·············Relacion: utc_adq (+/-) var_tiempo*nro_file= utc_pedestal time_Interval_p = 0.01 -n_perfiles_p = 100 +n_perfiles_p = 360 if utc_adq>utc_pedestal: nro_file = int((int(utc_adq) - int(utc_pedestal))/(time_Interval_p*n_perfiles_p)) ff_pedestal = list_pedestal[nro_file] @@ -129,7 +129,7 @@ list_adq = getfirstFilefromPath(path=dir_adq ,meta="D",ext=".hdf5") nro_file = nro_file #10 nro_key_perfil = nro_key_p -blocksPerFile = 100 +blocksPerFile = 360 wr_path = "/home/alex/Downloads/hdf5_wr/" # Lectura de archivos de adquisicion para adicion de azimuth for thisFile in range(len(list_adq)):