diff --git a/schainpy/__init__.py b/schainpy/__init__.py index d46a2d0..4654380 100644 --- a/schainpy/__init__.py +++ b/schainpy/__init__.py @@ -4,4 +4,4 @@ Created on Feb 7, 2012 @author $Author$ @version $Id$ ''' -__version__ = "2.3" +__version__ = '2.3' diff --git a/schainpy/model/data/BLTRheaderIO.py b/schainpy/model/data/BLTRheaderIO.py new file mode 100644 index 0000000..69e0d25 --- /dev/null +++ b/schainpy/model/data/BLTRheaderIO.py @@ -0,0 +1,404 @@ +''' + +$Author: murco $ +$Id: JROHeaderIO.py 151 2012-10-31 19:00:51Z murco $ +''' +import sys +import numpy +import copy +import datetime +from __builtin__ import None + +SPEED_OF_LIGHT = 299792458 +SPEED_OF_LIGHT = 3e8 + +FILE_STRUCTURE = numpy.dtype([ #HEADER 48bytes + ('FileMgcNumber',' vertical) + ('AntennaCoord',' endFp: + sys.stderr.write("Warning %s: Size value read from System Header is lower than it has to be\n" %fp.name) + return 0 + + if fp.tell() < endFp: + sys.stderr.write("Warning %s: Size value read from System Header size is greater than it has to be\n" %fp.name) + return 0 + + return 1 + + def write(self, fp): + + headerTuple = (self.RecMgcNumber, + self.RecCounter, + self.Off2StartNxtRec, + self.EpTimeStamp, + self.msCompTimeStamp, + self.ExpTagName, + self.ExpComment, + self.SiteLatDegrees, + self.SiteLongDegrees, + self.RTCgpsStatus, + self.TransmitFrec, + self.ReceiveFrec, + self.FirstOsciFrec, + self.Polarisation, + self.ReceiverFiltSett, + self.nModesInUse, + self.DualModeIndex, + self.DualModeRange, + self.nDigChannels, + self.SampResolution, + self.nRangeGatesSamp, + self.StartRangeSamp, + self.PRFhz, + self.Integrations, + self.nDataPointsTrsf, + self.nReceiveBeams, + self.nSpectAverages, + self.FFTwindowingInd, + self.BeamAngleAzim, + self.BeamAngleZen, + self.AntennaCoord, + self.RecPhaseCalibr, + self.RecAmpCalibr, + self.ReceiverGaindB) + +# self.size,self.nSamples, +# self.nProfiles, +# self.nChannels, +# self.adcResolution, +# self.pciDioBusWidth + + header = numpy.array(headerTuple,RECORD_STRUCTURE) + header.tofile(fp) + + return 1 + + +def get_dtype_index(numpy_dtype): + + index = None + + for i in range(len(NUMPY_DTYPE_LIST)): + if numpy_dtype == NUMPY_DTYPE_LIST[i]: + index = i + break + + return index + +def get_numpy_dtype(index): + + #dtype4 = numpy.dtype([('real','= self.nplots: break - + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) - + if showprofile: self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1) - + counter += 1 - + def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False, xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False): - + """ - + Input: dataOut : id : @@ -93,15 +93,15 @@ class CorrelationPlot(Figure): zmin : None, zmax : None """ - + if dataOut.flagNoData: return None - + if realtime: if not(isRealtime(utcdatatime = dataOut.utctime)): print 'Skipping this plot function' return - + if channelList == None: channelIndexList = dataOut.channelIndexList else: @@ -110,53 +110,53 @@ class CorrelationPlot(Figure): if channel not in dataOut.channelList: raise ValueError, "Channel %d is not in dataOut.channelList" channelIndexList.append(dataOut.channelList.index(channel)) - + factor = dataOut.normFactor lenfactor = factor.shape[1] x = dataOut.getLagTRange(1) y = dataOut.getHeiRange() - + z = copy.copy(dataOut.data_corr[:,:,0,:]) for i in range(dataOut.data_corr.shape[0]): - z[i,:,:] = z[i,:,:]/factor[i,:] + z[i,:,:] = z[i,:,:]/factor[i,:] zdB = numpy.abs(z) - + avg = numpy.average(z, axis=1) # avg = numpy.nanmean(z, axis=1) # noise = dataOut.noise/factor - + #thisDatetime = dataOut.datatime thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) - title = wintitle + " Correlation" + title = wintitle + " Correlation" xlabel = "Lag T (s)" ylabel = "Range (Km)" - + if not self.isConfig: - - nplots = dataOut.data_corr.shape[0] - + + nplots = dataOut.data_corr.shape[0] + self.setup(id=id, nplots=nplots, wintitle=wintitle, showprofile=showprofile, show=show) - + if xmin == None: xmin = numpy.nanmin(x) if xmax == None: xmax = numpy.nanmax(x) if ymin == None: ymin = numpy.nanmin(y) if ymax == None: ymax = numpy.nanmax(y) if zmin == None: zmin = 0 if zmax == None: zmax = 1 - + self.FTP_WEI = ftp_wei self.EXP_CODE = exp_code self.SUB_EXP_CODE = sub_exp_code self.PLOT_POS = plot_pos - + self.isConfig = True - + self.setWinTitle(title) - + for i in range(self.nplots): str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) title = "Channel %d and %d: : %s" %(dataOut.pairsList[i][0],dataOut.pairsList[i][1] , str_datetime) @@ -165,7 +165,7 @@ class CorrelationPlot(Figure): xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, xlabel=xlabel, ylabel=ylabel, title=title, ticksize=9, cblabel='') - + # if self.__showprofile: # axes = self.axesList[i*self.__nsubplots +1] # axes.pline(avgdB[i], y, @@ -173,15 +173,15 @@ class CorrelationPlot(Figure): # xlabel='dB', ylabel='', title='', # ytick_visible=False, # grid='x') -# +# # noiseline = numpy.repeat(noisedB[i], len(y)) # axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2) - + self.draw() - + self.save(figpath=figpath, figfile=figfile, save=save, ftp=ftp, wr_period=wr_period, - thisDatetime=thisDatetime) + thisDatetime=thisDatetime) diff --git a/schainpy/model/graphics/jroplot_data.py b/schainpy/model/graphics/jroplot_data.py index 8316adf..9e4caa7 100644 --- a/schainpy/model/graphics/jroplot_data.py +++ b/schainpy/model/graphics/jroplot_data.py @@ -15,10 +15,16 @@ from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator from schainpy.model.proc.jroproc_base import Operation from schainpy.utils import log -func = lambda x, pos: ('%s') %(datetime.datetime.fromtimestamp(x).strftime('%H:%M')) +jet_values = matplotlib.pyplot.get_cmap("jet", 100)(numpy.arange(100))[10:90] +blu_values = matplotlib.pyplot.get_cmap("seismic_r", 20)(numpy.arange(20))[10:15] +ncmap = matplotlib.colors.LinearSegmentedColormap.from_list("jro", numpy.vstack((blu_values, jet_values))) +matplotlib.pyplot.register_cmap(cmap=ncmap) -d1970 = datetime.datetime(1970, 1, 1) +func = lambda x, pos: '{}'.format(datetime.datetime.fromtimestamp(x).strftime('%H:%M')) +UT1970 = datetime.datetime(1970, 1, 1) - datetime.timedelta(seconds=time.timezone) + +CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'RdBu_r', 'seismic')] class PlotData(Operation, Process): ''' @@ -59,9 +65,7 @@ class PlotData(Operation, Process): self.zmin = kwargs.get('zmin', None) self.zmax = kwargs.get('zmax', None) self.zlimits = kwargs.get('zlimits', None) - self.xmin = kwargs.get('xmin', None) - if self.xmin is not None: - self.xmin += 5 + self.xmin = kwargs.get('xmin', None) self.xmax = kwargs.get('xmax', None) self.xrange = kwargs.get('xrange', 24) self.ymin = kwargs.get('ymin', None) @@ -83,6 +87,8 @@ class PlotData(Operation, Process): self.setup() + self.time_label = 'LT' if self.localtime else 'UTC' + if self.width is None: self.width = 8 @@ -106,6 +112,7 @@ class PlotData(Operation, Process): ax = fig.add_subplot(self.nrows, self.ncols, n+1) ax.tick_params(labelsize=8) ax.firsttime = True + ax.index = 0 self.axes.append(ax) if self.showprofile: cax = self.__add_axes(ax, size=size, pad=pad) @@ -121,6 +128,7 @@ class PlotData(Operation, Process): ax = fig.add_subplot(1, 1, 1) ax.tick_params(labelsize=8) ax.firsttime = True + ax.index = 0 self.figures.append(fig) self.axes.append(ax) if self.showprofile: @@ -136,6 +144,29 @@ class PlotData(Operation, Process): cmap.set_bad(self.bgcolor, 1.) self.cmaps.append(cmap) + for fig in self.figures: + fig.canvas.mpl_connect('key_press_event', self.event_key_press) + + def event_key_press(self, event): + ''' + ''' + + for ax in self.axes: + if ax == event.inaxes: + if event.key == 'down': + ax.index += 1 + elif event.key == 'up': + ax.index -= 1 + if ax.index < 0: + ax.index = len(CMAPS) - 1 + elif ax.index == len(CMAPS): + ax.index = 0 + cmap = CMAPS[ax.index] + ax.cbar.set_cmap(cmap) + ax.cbar.draw_all() + ax.plt.set_cmap(cmap) + ax.cbar.patch.figure.canvas.draw() + def __add_axes(self, ax, size='30%', pad='8%'): ''' Add new axes to the given figure @@ -145,6 +176,7 @@ class PlotData(Operation, Process): ax.figure.add_axes(nax) return nax + self.setup() def setup(self): ''' @@ -203,7 +235,7 @@ class PlotData(Operation, Process): if self.xaxis is 'time': dt = datetime.datetime.fromtimestamp(self.min_time) xmin = (datetime.datetime.combine(dt.date(), - datetime.time(int(self.xmin), 0, 0))-d1970).total_seconds() + datetime.time(int(self.xmin), 0, 0))-UT1970).total_seconds() else: xmin = self.xmin @@ -213,7 +245,7 @@ class PlotData(Operation, Process): if self.xaxis is 'time': dt = datetime.datetime.fromtimestamp(self.min_time) xmax = (datetime.datetime.combine(dt.date(), - datetime.time(int(self.xmax), 0, 0))-d1970).total_seconds() + datetime.time(int(self.xmax), 0, 0))-UT1970).total_seconds() else: xmax = self.xmax @@ -240,20 +272,20 @@ class PlotData(Operation, Process): self.pf_axes[n].grid(b=True, axis='x') [tick.set_visible(False) for tick in self.pf_axes[n].get_yticklabels()] if self.colorbar: - cb = plt.colorbar(ax.plt, ax=ax, pad=0.02) - cb.ax.tick_params(labelsize=8) + ax.cbar = plt.colorbar(ax.plt, ax=ax, pad=0.02, aspect=10) + ax.cbar.ax.tick_params(labelsize=8) if self.cb_label: - cb.set_label(self.cb_label, size=8) + ax.cbar.set_label(self.cb_label, size=8) elif self.cb_labels: - cb.set_label(self.cb_labels[n], size=8) - - ax.set_title('{} - {} UTC'.format( + ax.cbar.set_label(self.cb_labels[n], size=8) + + ax.set_title('{} - {} {}'.format( self.titles[n], - datetime.datetime.fromtimestamp(self.max_time).strftime('%H:%M:%S')), + datetime.datetime.fromtimestamp(self.max_time).strftime('%H:%M:%S'), + self.time_label), size=8) ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) - def __plot(self): ''' @@ -313,10 +345,15 @@ class PlotData(Operation, Process): while True: try: - self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK) + self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK) - self.min_time = self.data.times[0] - self.max_time = self.data.times[-1] + if self.localtime: + self.times = self.data.times - time.timezone + else: + self.times = self.data.times + + self.min_time = self.times[0] + self.max_time = self.times[-1] if self.isConfig is False: self.__setup() @@ -335,7 +372,6 @@ class PlotData(Operation, Process): if self.data: self.__plot() - class PlotSpectraData(PlotData): ''' Plot for Spectra data @@ -532,7 +568,7 @@ class PlotRTIData(PlotData): self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)] def plot(self): - self.x = self.data.times + self.x = self.times self.y = self.data.heights self.z = self.data[self.CODE] self.z = numpy.ma.masked_invalid(self.z) @@ -613,7 +649,7 @@ class PlotNoiseData(PlotData): def plot(self): - x = self.data.times + x = self.times xmin = self.min_time xmax = xmin+self.xrange*60*60 Y = self.data[self.CODE] @@ -681,7 +717,7 @@ class PlotSkyMapData(PlotData): def plot(self): - arrayParameters = numpy.concatenate([self.data['param'][t] for t in self.data.times]) + arrayParameters = numpy.concatenate([self.data['param'][t] for t in self.times]) error = arrayParameters[:,-1] indValid = numpy.where(error == 0)[0] finalMeteor = arrayParameters[indValid,:] @@ -727,6 +763,7 @@ class PlotParamData(PlotRTIData): self.nplots = self.nrows if self.showSNR: self.nrows += 1 + self.nplots += 1 self.ylabel = 'Height [Km]' self.titles = self.data.parameters \ @@ -736,7 +773,7 @@ class PlotParamData(PlotRTIData): def plot(self): self.data.normalize_heights() - self.x = self.data.times + self.x = self.times self.y = self.data.heights if self.showSNR: self.z = numpy.concatenate( @@ -779,4 +816,4 @@ class PlotOuputData(PlotParamData): ''' CODE = 'output' - colormap = 'seismic' \ No newline at end of file + colormap = 'seismic' diff --git a/schainpy/model/graphics/jroplot_heispectra.py b/schainpy/model/graphics/jroplot_heispectra.py index b0c9c5b..f8d4512 100644 --- a/schainpy/model/graphics/jroplot_heispectra.py +++ b/schainpy/model/graphics/jroplot_heispectra.py @@ -11,80 +11,79 @@ from figure import Figure, isRealtime from plotting_codes import * class SpectraHeisScope(Figure): - - + + isConfig = None __nsubplots = None - + WIDTHPROF = None HEIGHTPROF = None PREFIX = 'spc' - - def __init__(self, **kwargs): - - Figure.__init__(self, **kwargs) + + def __init__(self): + self.isConfig = False self.__nsubplots = 1 - + self.WIDTH = 230 self.HEIGHT = 250 self.WIDTHPROF = 120 self.HEIGHTPROF = 0 self.counter_imagwr = 0 - + self.PLOT_CODE = SPEC_CODE - + def getSubplots(self): - + ncol = int(numpy.sqrt(self.nplots)+0.9) nrow = int(self.nplots*1./ncol + 0.9) - + return nrow, ncol - + def setup(self, id, nplots, wintitle, show): - + showprofile = False self.__showprofile = showprofile self.nplots = nplots - + ncolspan = 1 colspan = 1 if showprofile: ncolspan = 3 colspan = 2 self.__nsubplots = 2 - + self.createFigure(id = id, wintitle = wintitle, widthplot = self.WIDTH + self.WIDTHPROF, heightplot = self.HEIGHT + self.HEIGHTPROF, show = show) - + nrow, ncol = self.getSubplots() - + counter = 0 for y in range(nrow): for x in range(ncol): - + if counter >= self.nplots: break - + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) - + if showprofile: self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1) - + counter += 1 - + def run(self, dataOut, id, wintitle="", channelList=None, xmin=None, xmax=None, ymin=None, ymax=None, save=False, figpath='./', figfile=None, ftp=False, wr_period=1, show=True, server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): - + """ - + Input: dataOut : id : @@ -95,12 +94,12 @@ class SpectraHeisScope(Figure): ymin : None, ymax : None, """ - + if dataOut.realtime: if not(isRealtime(utcdatatime = dataOut.utctime)): print 'Skipping this plot function' return - + if channelList == None: channelIndexList = dataOut.channelIndexList else: @@ -109,9 +108,9 @@ class SpectraHeisScope(Figure): if channel not in dataOut.channelList: raise ValueError, "Channel %d is not in dataOut.channelList" channelIndexList.append(dataOut.channelList.index(channel)) - + # x = dataOut.heightList - c = 3E8 + c = 3E8 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] #deberia cambiar para el caso de 1Mhz y 100KHz x = numpy.arange(-1*dataOut.nHeights/2.,dataOut.nHeights/2.)*(c/(2*deltaHeight*dataOut.nHeights*1000)) @@ -123,7 +122,7 @@ class SpectraHeisScope(Figure): data = dataOut.data_spc / factor datadB = 10.*numpy.log10(data) y = datadB - + #thisDatetime = dataOut.datatime thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) @@ -131,29 +130,29 @@ class SpectraHeisScope(Figure): #para 1Mhz descomentar la siguiente linea #xlabel = "Frequency x 10000" ylabel = "Intensity (dB)" - + if not self.isConfig: nplots = len(channelIndexList) - + self.setup(id=id, nplots=nplots, wintitle=wintitle, show=show) - + if xmin == None: xmin = numpy.nanmin(x) if xmax == None: xmax = numpy.nanmax(x) if ymin == None: ymin = numpy.nanmin(y) if ymax == None: ymax = numpy.nanmax(y) - + self.FTP_WEI = ftp_wei self.EXP_CODE = exp_code self.SUB_EXP_CODE = sub_exp_code self.PLOT_POS = plot_pos - + self.isConfig = True - + self.setWinTitle(title) - + for i in range(len(self.axesList)): ychannel = y[i,:] str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) @@ -162,10 +161,10 @@ class SpectraHeisScope(Figure): axes.pline(x, ychannel, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, xlabel=xlabel, ylabel=ylabel, title=title, grid='both') - - + + self.draw() - + self.save(figpath=figpath, figfile=figfile, save=save, @@ -174,18 +173,18 @@ class SpectraHeisScope(Figure): thisDatetime=thisDatetime) class RTIfromSpectraHeis(Figure): - + isConfig = None __nsubplots = None PREFIX = 'rtinoise' - - def __init__(self, **kwargs): - Figure.__init__(self, **kwargs) + + def __init__(self): + self.timerange = 24*60*60 self.isConfig = False self.__nsubplots = 1 - + self.WIDTH = 820 self.HEIGHT = 200 self.WIDTHPROF = 120 @@ -194,43 +193,43 @@ class RTIfromSpectraHeis(Figure): self.xdata = None self.ydata = None self.figfile = None - + self.PLOT_CODE = RTI_CODE - + 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 run(self, dataOut, id, wintitle="", channelList=None, showprofile='True', xmin=None, xmax=None, ymin=None, ymax=None, timerange=None, save=False, figpath='./', figfile=None, ftp=False, wr_period=1, show=True, server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): - + if channelList == None: channelIndexList = dataOut.channelIndexList channelList = dataOut.channelList @@ -240,86 +239,86 @@ class RTIfromSpectraHeis(Figure): if channel not in dataOut.channelList: raise ValueError, "Channel %d is not in dataOut.channelList" channelIndexList.append(dataOut.channelList.index(channel)) - + if timerange != None: self.timerange = timerange - + x = dataOut.getTimeRange() y = dataOut.getHeiRange() - + factor = dataOut.normFactor data = dataOut.data_spc / factor data = numpy.average(data,axis=1) datadB = 10*numpy.log10(data) - + # factor = dataOut.normFactor # noise = dataOut.getNoise()/factor # noisedB = 10*numpy.log10(noise) - + #thisDatetime = dataOut.datatime thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + " RTI: %s" %(thisDatetime.strftime("%d-%b-%Y")) xlabel = "Local Time" ylabel = "Intensity (dB)" - + if not self.isConfig: - + nplots = 1 - + self.setup(id=id, nplots=nplots, wintitle=wintitle, showprofile=showprofile, show=show) - + self.tmin, self.tmax = self.getTimeLim(x, xmin, xmax) - + if ymin == None: ymin = numpy.nanmin(datadB) if ymax == None: ymax = numpy.nanmax(datadB) - + self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") self.isConfig = True self.figfile = figfile self.xdata = numpy.array([]) self.ydata = numpy.array([]) - + self.FTP_WEI = ftp_wei self.EXP_CODE = exp_code self.SUB_EXP_CODE = sub_exp_code self.PLOT_POS = plot_pos - + self.setWinTitle(title) - - + + # title = "RTI %s" %(thisDatetime.strftime("%d-%b-%Y")) title = "RTI - %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) - + legendlabels = ["channel %d"%idchannel for idchannel in channelList] axes = self.axesList[0] - + self.xdata = numpy.hstack((self.xdata, x[0:1])) - + if len(self.ydata)==0: self.ydata = datadB[channelIndexList].reshape(-1,1) else: self.ydata = numpy.hstack((self.ydata, datadB[channelIndexList].reshape(-1,1))) - - + + axes.pmultilineyaxis(x=self.xdata, y=self.ydata, xmin=self.tmin, xmax=self.tmax, ymin=ymin, ymax=ymax, xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='.', markersize=8, linestyle="solid", grid='both', XAxisAsTime=True ) - + self.draw() - + update_figfile = False - + if dataOut.ltctime >= self.tmax: self.counter_imagwr = wr_period self.isConfig = False update_figfile = True - + self.save(figpath=figpath, figfile=figfile, save=save, diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index ffc1924..4e810c6 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -6,6 +6,217 @@ from figure import Figure, isRealtime, isTimeInHourRange from plotting_codes import * +class FitGauPlot(Figure): + + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'fitgau' + + def __init__(self, **kwargs): + Figure.__init__(self, **kwargs) + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 250 + self.HEIGHT = 250 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + self.PLOT_CODE = SPEC_CODE + + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + + self.__xfilter_ena = False + self.__yfilter_ena = False + + def getSubplots(self): + + ncol = int(numpy.sqrt(self.nplots)+0.9) + nrow = int(self.nplots*1./ncol + 0.9) + + return nrow, ncol + + def setup(self, id, nplots, wintitle, showprofile=True, show=True): + + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + if showprofile: + ncolspan = 3 + colspan = 2 + self.__nsubplots = 2 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH + self.WIDTHPROF, + heightplot = self.HEIGHT + self.HEIGHTPROF, + show=show) + + nrow, ncol = self.getSubplots() + + counter = 0 + for y in range(nrow): + for x in range(ncol): + + if counter >= self.nplots: + break + + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) + + if showprofile: + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1) + + counter += 1 + + def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True, + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False, + xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 1): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : None + """ + if realtime: + if not(isRealtime(utcdatatime = dataOut.utctime)): + print 'Skipping this plot function' + return + + if channelList == None: + channelIndexList = dataOut.channelIndexList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" %channel + channelIndexList.append(dataOut.channelList.index(channel)) + +# if normFactor is None: +# factor = dataOut.normFactor +# else: +# factor = normFactor + if xaxis == "frequency": + x = dataOut.spc_range[0] + xlabel = "Frequency (kHz)" + + elif xaxis == "time": + x = dataOut.spc_range[1] + xlabel = "Time (ms)" + + else: + x = dataOut.spc_range[2] + xlabel = "Velocity (m/s)" + + ylabel = "Range (Km)" + + y = dataOut.getHeiRange() + + z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor + print 'GausSPC', z[0,32,10:40] + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + zdB = 10*numpy.log10(z) + + avg = numpy.average(z, axis=1) + avgdB = 10*numpy.log10(avg) + + noise = dataOut.spc_noise + noisedB = 10*numpy.log10(noise) + + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) + title = wintitle + " Spectra" + if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)): + title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith) + + if not self.isConfig: + + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + if zmin == None: zmin = numpy.floor(numpy.nanmin(noisedB)) - 3 + if zmax == None: zmax = numpy.ceil(numpy.nanmax(avgdB)) + 3 + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(self.nplots): + index = channelIndexList[i] + str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) + title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[index], noisedB[index], str_datetime) + if len(dataOut.beam.codeList) != 0: + title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[index], noisedB[index], dataOut.beam.azimuthList[index], dataOut.beam.zenithList[index], str_datetime) + + axes = self.axesList[i*self.__nsubplots] + axes.pcolor(x, y, zdB[index,:,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, colormap=colormap, + ticksize=9, cblabel='') + + if self.__showprofile: + axes = self.axesList[i*self.__nsubplots +1] + axes.pline(avgdB[index,:], y, + xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax, + xlabel='dB', ylabel='', title='', + ytick_visible=False, + grid='x') + + noiseline = numpy.repeat(noisedB[index], len(y)) + axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2) + + self.draw() + + if figfile == None: + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + name = str_datetime + if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)): + name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith) + figfile = self.getFilename(name) + + self.save(figpath=figpath, + figfile=figfile, + save=save, + ftp=ftp, + wr_period=wr_period, + thisDatetime=thisDatetime) + + + class MomentsPlot(Figure): isConfig = None @@ -445,10 +656,9 @@ class WindProfilerPlot(Figure): # tmin = None # tmax = None - - x = dataOut.getTimeRange1(dataOut.outputInterval) - y = dataOut.heightList - z = dataOut.data_output.copy() + x = dataOut.getTimeRange1(dataOut.paramInterval) + y = dataOut.heightList + z = dataOut.data_output.copy() nplots = z.shape[0] #Number of wind dimensions estimated nplotsw = nplots @@ -558,7 +768,7 @@ class WindProfilerPlot(Figure): thisDatetime=thisDatetime, update_figfile=update_figfile) - if dataOut.ltctime + dataOut.outputInterval >= self.xmax: + if dataOut.ltctime + dataOut.paramInterval >= self.xmax: self.counter_imagwr = wr_period self.isConfig = False update_figfile = True @@ -635,12 +845,12 @@ class ParametersPlot(Figure): counter += 1 - def run(self, dataOut, id, wintitle="", channelList=None, paramIndex = 0, colormap=True, + def run(self, dataOut, id, wintitle="", channelList=None, paramIndex = 0, colormap="jet", xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, timerange=None, showSNR=False, SNRthresh = -numpy.inf, SNRmin=None, SNRmax=None, save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True, server=None, folder=None, username=None, password=None, - ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, HEIGHT=None): """ Input: @@ -656,12 +866,11 @@ class ParametersPlot(Figure): zmin : None, zmax : None """ - - if colormap: - colormap="jet" - else: - colormap="RdBu_r" - + + if HEIGHT is not None: + self.HEIGHT = HEIGHT + + if not isTimeInHourRange(dataOut.datatime, xmin, xmax): return @@ -906,25 +1115,23 @@ class Parameters1Plot(Figure): x = dataOut.getTimeRange1(dataOut.paramInterval) y = dataOut.heightList - z = data_param[channelIndexList,parameterIndex,:].copy() - zRange = dataOut.abscissaList -# nChannels = z.shape[0] #Number of wind dimensions estimated -# thisDatetime = dataOut.datatime + if dataOut.data_param.ndim == 3: + z = dataOut.data_param[channelIndexList,parameterIndex,:] + else: + z = dataOut.data_param[channelIndexList,:] if dataOut.data_SNR is not None: - SNRarray = dataOut.data_SNR[channelIndexList,:] - SNRdB = 10*numpy.log10(SNRarray) -# SNRavgdB = 10*numpy.log10(SNRavg) - ind = numpy.where(SNRdB < 10**(SNRthresh/10)) - z[ind] = numpy.nan + if dataOut.data_SNR.ndim == 2: + SNRavg = numpy.average(dataOut.data_SNR, axis=0) + else: + SNRavg = dataOut.data_SNR + SNRdB = 10*numpy.log10(SNRavg) thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y")) xlabel = "" - ylabel = "Range (Km)" - - if (SNR and not onlySNR): nplots = 2*nplots + ylabel = "Range (Km)" if onlyPositive: colormap = "jet" @@ -943,8 +1150,8 @@ class Parameters1Plot(Figure): if ymin == None: ymin = numpy.nanmin(y) if ymax == None: ymax = numpy.nanmax(y) - if zmin == None: zmin = numpy.nanmin(zRange) - if zmax == None: zmax = numpy.nanmax(zRange) + if zmin == None: zmin = numpy.nanmin(z) + if zmax == None: zmax = numpy.nanmax(z) if SNR: if SNRmin == None: SNRmin = numpy.nanmin(SNRdB) @@ -994,19 +1201,18 @@ class Parameters1Plot(Figure): xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap, ticksize=9, cblabel=zlabel, cbsize="1%") - if SNR: - title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) - axes = self.axesList[(j)*self.__nsubplots] - if not onlySNR: - axes = self.axesList[(j + 1)*self.__nsubplots] - - axes = self.axesList[(j + nGraphsByChannel-1)] + if SNR: + title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[(j)*self.__nsubplots] + if not onlySNR: + axes = self.axesList[(j + 1)*self.__nsubplots] - z1 = SNRdB[i,:].reshape((1,-1)) - axes.pcolorbuffer(x, y, z1, - xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, - xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap="jet", - ticksize=9, cblabel=zlabel, cbsize="1%") + axes = self.axesList[(j + nGraphsByChannel-1)] + z1 = SNRdB.reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap="jet", + ticksize=9, cblabel=zlabel, cbsize="1%") diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 87efee7..eafac9f 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -87,7 +87,7 @@ class SpectraPlot(Figure): save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False, - xaxis="velocity", **kwargs): + xaxis="frequency", colormap='jet', normFactor=None): """ @@ -104,9 +104,6 @@ class SpectraPlot(Figure): zmin : None, zmax : None """ - - colormap = kwargs.get('colormap','jet') - if realtime: if not(isRealtime(utcdatatime = dataOut.utctime)): print 'Skipping this plot function' @@ -121,8 +118,10 @@ class SpectraPlot(Figure): raise ValueError, "Channel %d is not in dataOut.channelList" %channel channelIndexList.append(dataOut.channelList.index(channel)) - factor = dataOut.normFactor - + if normFactor is None: + factor = dataOut.normFactor + else: + factor = normFactor if xaxis == "frequency": x = dataOut.getFreqRange(1)/1000. xlabel = "Frequency (kHz)" @@ -283,7 +282,7 @@ class CrossSpectraPlot(Figure): save=False, figpath='./', figfile=None, ftp=False, wr_period=1, power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True, server=None, folder=None, username=None, password=None, - ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, xaxis='frequency'): """ @@ -316,8 +315,11 @@ class CrossSpectraPlot(Figure): if len(pairsIndexList) > 4: pairsIndexList = pairsIndexList[0:4] - - factor = dataOut.normFactor + + if normFactor is None: + factor = dataOut.normFactor + else: + factor = normFactor x = dataOut.getVelRange(1) y = dataOut.getHeiRange() z = dataOut.data_spc[:,:,:]/factor @@ -518,10 +520,10 @@ class RTIPlot(Figure): def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True', xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, - timerange=None, + timerange=None, colormap='jet', save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True, server=None, folder=None, username=None, password=None, - ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, **kwargs): + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, normFactor=None, HEIGHT=None): """ @@ -539,7 +541,10 @@ class RTIPlot(Figure): zmax : None """ - colormap = kwargs.get('colormap', 'jet') + #colormap = kwargs.get('colormap', 'jet') + if HEIGHT is not None: + self.HEIGHT = HEIGHT + if not isTimeInHourRange(dataOut.datatime, xmin, xmax): return @@ -552,20 +557,21 @@ class RTIPlot(Figure): raise ValueError, "Channel %d is not in dataOut.channelList" channelIndexList.append(dataOut.channelList.index(channel)) - if hasattr(dataOut, 'normFactor'): + if normFactor is None: factor = dataOut.normFactor else: - factor = 1 + factor = normFactor # factor = dataOut.normFactor x = dataOut.getTimeRange() y = dataOut.getHeiRange() -# z = dataOut.data_spc/factor -# z = numpy.where(numpy.isfinite(z), z, numpy.NAN) -# avg = numpy.average(z, axis=1) -# avgdB = 10.*numpy.log10(avg) - avgdB = dataOut.getPower() + z = dataOut.data_spc/factor + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + avg = numpy.average(z, axis=1) + avgdB = 10.*numpy.log10(avg) + # avgdB = dataOut.getPower() + thisDatetime = dataOut.datatime # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) @@ -1112,6 +1118,7 @@ class Noise(Figure): PREFIX = 'noise' + def __init__(self, **kwargs): Figure.__init__(self, **kwargs) self.timerange = 24*60*60 diff --git a/schainpy/model/graphics/mpldriver.py b/schainpy/model/graphics/mpldriver.py index b84fb1c..0b3b227 100644 --- a/schainpy/model/graphics/mpldriver.py +++ b/schainpy/model/graphics/mpldriver.py @@ -4,7 +4,7 @@ import sys import matplotlib if 'linux' in sys.platform: - matplotlib.use("TKAgg") + matplotlib.use("GTK3Agg") if 'darwin' in sys.platform: matplotlib.use('TKAgg') @@ -88,6 +88,8 @@ def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False): rowspan=rowspan, polar=polar) + axes.grid(True) + matplotlib.pyplot.ion() return axes @@ -174,7 +176,7 @@ def set_linedata(ax, x, y, idline): def pline(iplot, x, y, xlabel='', ylabel='', title=''): - ax = iplot.axes + ax = iplot.get_axes() printLabels(ax, xlabel, ylabel, title) @@ -204,7 +206,7 @@ def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax, z = numpy.ma.masked_invalid(z) cmap=matplotlib.pyplot.get_cmap(colormap) - cmap.set_bad('black', 1.) + cmap.set_bad('white',1.) imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=cmap) cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb) cb.set_label(cblabel) @@ -239,21 +241,32 @@ def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax, ax.xaxis.set_major_formatter(FuncFormatter(func)) ax.xaxis.set_major_locator(LinearLocator(7)) + ax.grid(True) matplotlib.pyplot.ion() return imesh def pcolor(imesh, z, xlabel='', ylabel='', title=''): + z = numpy.ma.masked_invalid(z) + + cmap=matplotlib.pyplot.get_cmap('jet') + cmap.set_bad('white',1.) + z = z.T - ax = imesh.axes + ax = imesh.get_axes() printLabels(ax, xlabel, ylabel, title) imesh.set_array(z.ravel()) + ax.grid(True) + def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'): printLabels(ax, xlabel, ylabel, title) - + z = numpy.ma.masked_invalid(z) + cmap=matplotlib.pyplot.get_cmap(colormap) + cmap.set_bad('white',1.) ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap)) + ax.grid(True) def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'): @@ -262,12 +275,13 @@ def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', col ax.collections.remove(ax.collections[0]) z = numpy.ma.masked_invalid(z) - + cmap=matplotlib.pyplot.get_cmap(colormap) - cmap.set_bad('black', 1.) - + cmap.set_bad('white',1.) ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=cmap) + ax.grid(True) + def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None, ticksize=9, xtick_visible=True, ytick_visible=True, @@ -326,7 +340,7 @@ def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', tit def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''): - ax = iplot.axes + ax = iplot.get_axes() printLabels(ax, xlabel, ylabel, title) @@ -403,8 +417,7 @@ def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='' def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''): - ax = iplot.axes - + ax = iplot.get_axes() printLabels(ax, xlabel, ylabel, title) for i in range(len(ax.lines)): @@ -425,7 +438,7 @@ def createPolar(ax, x, y, # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11') # ax.text(0, 50, ylabel, rotation='vertical', va ='center', ha = 'left' ,size='11') # ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical') - ax.yaxis.labelpad = 40 + ax.yaxis.labelpad = 230 printLabels(ax, xlabel, ylabel, title) iplot = ax.lines[-1] @@ -449,7 +462,7 @@ def createPolar(ax, x, y, def polar(iplot, x, y, xlabel='', ylabel='', title=''): - ax = iplot.axes + ax = iplot.get_axes() # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11') printLabels(ax, xlabel, ylabel, title) diff --git a/schainpy/model/graphics/mpldriver2.py b/schainpy/model/graphics/mpldriver2.py new file mode 100644 index 0000000..d6cbfd1 --- /dev/null +++ b/schainpy/model/graphics/mpldriver2.py @@ -0,0 +1,469 @@ +import numpy +import datetime +import sys +import matplotlib + +if 'linux' in sys.platform: + matplotlib.use("TKAgg") + +if 'darwin' in sys.platform: + matplotlib.use('TKAgg') +#Qt4Agg', 'GTK', 'GTKAgg', 'ps', 'agg', 'cairo', 'MacOSX', 'GTKCairo', 'WXAgg', 'template', 'TkAgg', 'GTK3Cairo', 'GTK3Agg', 'svg', 'WebAgg', 'CocoaAgg', 'emf', 'gdk', 'WX' +import matplotlib.pyplot + +from mpl_toolkits.axes_grid1 import make_axes_locatable +from matplotlib.ticker import FuncFormatter, LinearLocator + +########################################### +#Actualizacion de las funciones del driver +########################################### + +jet_values = matplotlib.pyplot.get_cmap("jet", 100)(numpy.arange(100))[10:90] +blu_values = matplotlib.pyplot.get_cmap("seismic_r", 20)(numpy.arange(20))[10:15] +ncmap = matplotlib.colors.LinearSegmentedColormap.from_list("jro", numpy.vstack((blu_values, jet_values))) +matplotlib.pyplot.register_cmap(cmap=ncmap) + +def createFigure(id, wintitle, width, height, facecolor="w", show=True, dpi = 80): + + matplotlib.pyplot.ioff() + + fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor, figsize=(1.0*width/dpi, 1.0*height/dpi)) + fig.canvas.manager.set_window_title(wintitle) +# fig.canvas.manager.resize(width, height) + matplotlib.pyplot.ion() + + + if show: + matplotlib.pyplot.show() + + return fig + +def closeFigure(show=False, fig=None): + +# matplotlib.pyplot.ioff() +# matplotlib.pyplot.pause(0) + + if show: + matplotlib.pyplot.show() + + if fig != None: + matplotlib.pyplot.close(fig) +# matplotlib.pyplot.pause(0) +# matplotlib.pyplot.ion() + + return + + matplotlib.pyplot.close("all") +# matplotlib.pyplot.pause(0) +# matplotlib.pyplot.ion() + + return + +def saveFigure(fig, filename): + +# matplotlib.pyplot.ioff() + fig.savefig(filename, dpi=matplotlib.pyplot.gcf().dpi) +# matplotlib.pyplot.ion() + +def clearFigure(fig): + + fig.clf() + +def setWinTitle(fig, title): + + fig.canvas.manager.set_window_title(title) + +def setTitle(fig, title): + + fig.suptitle(title) + +def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False): + + matplotlib.pyplot.ioff() + matplotlib.pyplot.figure(fig.number) + axes = matplotlib.pyplot.subplot2grid((nrow, ncol), + (xpos, ypos), + colspan=colspan, + rowspan=rowspan, + polar=polar) + + axes.grid(True) + matplotlib.pyplot.ion() + return axes + +def setAxesText(ax, text): + + ax.annotate(text, + xy = (.1, .99), + xycoords = 'figure fraction', + horizontalalignment = 'left', + verticalalignment = 'top', + fontsize = 10) + +def printLabels(ax, xlabel, ylabel, title): + + ax.set_xlabel(xlabel, size=11) + ax.set_ylabel(ylabel, size=11) + ax.set_title(title, size=8) + +def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', + ticksize=9, xtick_visible=True, ytick_visible=True, + nxticks=4, nyticks=10, + grid=None,color='blue'): + + """ + + Input: + grid : None, 'both', 'x', 'y' + """ + + matplotlib.pyplot.ioff() + + ax.set_xlim([xmin,xmax]) + ax.set_ylim([ymin,ymax]) + + printLabels(ax, xlabel, ylabel, title) + + ###################################################### + if (xmax-xmin)<=1: + xtickspos = numpy.linspace(xmin,xmax,nxticks) + xtickspos = numpy.array([float("%.1f"%i) for i in xtickspos]) + ax.set_xticks(xtickspos) + else: + xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin) +# xtickspos = numpy.arange(nxticks)*float(xmax-xmin)/float(nxticks) + int(xmin) + ax.set_xticks(xtickspos) + + for tick in ax.get_xticklabels(): + tick.set_visible(xtick_visible) + + for tick in ax.xaxis.get_major_ticks(): + tick.label.set_fontsize(ticksize) + + ###################################################### + for tick in ax.get_yticklabels(): + tick.set_visible(ytick_visible) + + for tick in ax.yaxis.get_major_ticks(): + tick.label.set_fontsize(ticksize) + + ax.plot(x, y, color=color) + iplot = ax.lines[-1] + + ###################################################### + if '0.' in matplotlib.__version__[0:2]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + + if '1.0.' in matplotlib.__version__[0:4]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + + if grid != None: + ax.grid(b=True, which='major', axis=grid) + + matplotlib.pyplot.tight_layout() + + matplotlib.pyplot.ion() + + return iplot + +def set_linedata(ax, x, y, idline): + + ax.lines[idline].set_data(x,y) + +def pline(iplot, x, y, xlabel='', ylabel='', title=''): + + ax = iplot.get_axes() + + printLabels(ax, xlabel, ylabel, title) + + set_linedata(ax, x, y, idline=0) + +def addpline(ax, x, y, color, linestyle, lw): + + ax.plot(x,y,color=color,linestyle=linestyle,lw=lw) + + +def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax, + xlabel='', ylabel='', title='', ticksize = 9, + colormap='jet',cblabel='', cbsize="5%", + XAxisAsTime=False): + + matplotlib.pyplot.ioff() + + divider = make_axes_locatable(ax) + ax_cb = divider.new_horizontal(size=cbsize, pad=0.05) + fig = ax.get_figure() + fig.add_axes(ax_cb) + + ax.set_xlim([xmin,xmax]) + ax.set_ylim([ymin,ymax]) + + printLabels(ax, xlabel, ylabel, title) + + z = numpy.ma.masked_invalid(z) + cmap=matplotlib.pyplot.get_cmap(colormap) + cmap.set_bad('white',1.) + imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=cmap) + cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb) + cb.set_label(cblabel) + +# for tl in ax_cb.get_yticklabels(): +# tl.set_visible(True) + + for tick in ax.yaxis.get_major_ticks(): + tick.label.set_fontsize(ticksize) + + for tick in ax.xaxis.get_major_ticks(): + tick.label.set_fontsize(ticksize) + + for tick in cb.ax.get_yticklabels(): + tick.set_fontsize(ticksize) + + ax_cb.yaxis.tick_right() + + if '0.' in matplotlib.__version__[0:2]: + print "The matplotlib version has to be updated to 1.1 or newer" + return imesh + + if '1.0.' in matplotlib.__version__[0:4]: + print "The matplotlib version has to be updated to 1.1 or newer" + return imesh + + matplotlib.pyplot.tight_layout() + + if XAxisAsTime: + + func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S")) + ax.xaxis.set_major_formatter(FuncFormatter(func)) + ax.xaxis.set_major_locator(LinearLocator(7)) + ax.grid(True) + matplotlib.pyplot.ion() + return imesh + +def pcolor(imesh, z, xlabel='', ylabel='', title=''): + + z = z.T + ax = imesh.get_axes() + printLabels(ax, xlabel, ylabel, title) + imesh.set_array(z.ravel()) + ax.grid(True) + +def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'): + + printLabels(ax, xlabel, ylabel, title) + ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap)) + ax.grid(True) + +def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'): + + printLabels(ax, xlabel, ylabel, title) + + ax.collections.remove(ax.collections[0]) + + z = numpy.ma.masked_invalid(z) + + cmap=matplotlib.pyplot.get_cmap(colormap) + cmap.set_bad('white',1.) + + ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=cmap) + ax.grid(True) + +def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None, + ticksize=9, xtick_visible=True, ytick_visible=True, + nxticks=4, nyticks=10, + grid=None): + + """ + + Input: + grid : None, 'both', 'x', 'y' + """ + + matplotlib.pyplot.ioff() + + lines = ax.plot(x.T, y) + leg = ax.legend(lines, legendlabels, loc='upper right') + leg.get_frame().set_alpha(0.5) + ax.set_xlim([xmin,xmax]) + ax.set_ylim([ymin,ymax]) + printLabels(ax, xlabel, ylabel, title) + + xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin) + ax.set_xticks(xtickspos) + + for tick in ax.get_xticklabels(): + tick.set_visible(xtick_visible) + + for tick in ax.xaxis.get_major_ticks(): + tick.label.set_fontsize(ticksize) + + for tick in ax.get_yticklabels(): + tick.set_visible(ytick_visible) + + for tick in ax.yaxis.get_major_ticks(): + tick.label.set_fontsize(ticksize) + + iplot = ax.lines[-1] + + if '0.' in matplotlib.__version__[0:2]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + + if '1.0.' in matplotlib.__version__[0:4]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + + if grid != None: + ax.grid(b=True, which='major', axis=grid) + + matplotlib.pyplot.tight_layout() + + matplotlib.pyplot.ion() + + return iplot + + +def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''): + + ax = iplot.get_axes() + + printLabels(ax, xlabel, ylabel, title) + + for i in range(len(ax.lines)): + line = ax.lines[i] + line.set_data(x[i,:],y) + +def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None, + ticksize=9, xtick_visible=True, ytick_visible=True, + nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None", + grid=None, XAxisAsTime=False): + + """ + + Input: + grid : None, 'both', 'x', 'y' + """ + + matplotlib.pyplot.ioff() + +# lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle) + lines = ax.plot(x, y.T) +# leg = ax.legend(lines, legendlabels, loc=2, bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \ +# handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.) + + leg = ax.legend(lines, legendlabels, + loc='upper right', bbox_to_anchor=(1.16, 1), borderaxespad=0) + + for label in leg.get_texts(): label.set_fontsize(9) + + ax.set_xlim([xmin,xmax]) + ax.set_ylim([ymin,ymax]) + printLabels(ax, xlabel, ylabel, title) + +# xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin) +# ax.set_xticks(xtickspos) + + for tick in ax.get_xticklabels(): + tick.set_visible(xtick_visible) + + for tick in ax.xaxis.get_major_ticks(): + tick.label.set_fontsize(ticksize) + + for tick in ax.get_yticklabels(): + tick.set_visible(ytick_visible) + + for tick in ax.yaxis.get_major_ticks(): + tick.label.set_fontsize(ticksize) + + iplot = ax.lines[-1] + + if '0.' in matplotlib.__version__[0:2]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + + if '1.0.' in matplotlib.__version__[0:4]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + + if grid != None: + ax.grid(b=True, which='major', axis=grid) + + matplotlib.pyplot.tight_layout() + + if XAxisAsTime: + + func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S")) + ax.xaxis.set_major_formatter(FuncFormatter(func)) + ax.xaxis.set_major_locator(LinearLocator(7)) + + matplotlib.pyplot.ion() + + return iplot + +def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''): + + ax = iplot.get_axes() + + printLabels(ax, xlabel, ylabel, title) + + for i in range(len(ax.lines)): + line = ax.lines[i] + line.set_data(x,y[i,:]) + +def createPolar(ax, x, y, + xlabel='', ylabel='', title='', ticksize = 9, + colormap='jet',cblabel='', cbsize="5%", + XAxisAsTime=False): + + matplotlib.pyplot.ioff() + + ax.plot(x,y,'bo', markersize=5) +# ax.set_rmax(90) + ax.set_ylim(0,90) + ax.set_yticks(numpy.arange(0,90,20)) +# ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11') +# ax.text(0, 50, ylabel, rotation='vertical', va ='center', ha = 'left' ,size='11') +# ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical') + ax.yaxis.labelpad = 230 + printLabels(ax, xlabel, ylabel, title) + iplot = ax.lines[-1] + + if '0.' in matplotlib.__version__[0:2]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + + if '1.0.' in matplotlib.__version__[0:4]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + +# if grid != None: +# ax.grid(b=True, which='major', axis=grid) + + matplotlib.pyplot.tight_layout() + + matplotlib.pyplot.ion() + + + return iplot + +def polar(iplot, x, y, xlabel='', ylabel='', title=''): + + ax = iplot.get_axes() + +# ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11') + printLabels(ax, xlabel, ylabel, title) + + set_linedata(ax, x, y, idline=0) + +def draw(fig): + + if type(fig) == 'int': + raise ValueError, "Error drawing: Fig parameter should be a matplotlib figure object figure" + + fig.canvas.draw() + +def pause(interval=0.000001): + + matplotlib.pyplot.pause(interval) diff --git a/schainpy/model/io/MIRAtest.py b/schainpy/model/io/MIRAtest.py new file mode 100644 index 0000000..c48bc8b --- /dev/null +++ b/schainpy/model/io/MIRAtest.py @@ -0,0 +1,321 @@ +import os, sys +import glob +import fnmatch +import datetime +import time +import re +import h5py +import numpy +import matplotlib.pyplot as plt + +import pylab as plb +from scipy.optimize import curve_fit +from scipy import asarray as ar,exp +from scipy import stats + +from duplicity.path import Path +from numpy.ma.core import getdata + +SPEED_OF_LIGHT = 299792458 +SPEED_OF_LIGHT = 3e8 + +try: + from gevent import sleep +except: + from time import sleep + +from schainpy.model.data.jrodata import Spectra +#from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader +from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation +#from schainpy.model.io.jroIO_bltr import BLTRReader +from numpy import imag, shape, NaN + + +startFp = open('/home/erick/Documents/MIRA35C/20160117/20160117_0000.zspc',"rb") + + +FILE_HEADER = numpy.dtype([ #HEADER 1024bytes + ('Hname',numpy.str_,32), #Original file name + ('Htime',numpy.str_,32), #Date and time when the file was created + ('Hoper',numpy.str_,64), #Name of operator who created the file + ('Hplace',numpy.str_,128), #Place where the measurements was carried out + ('Hdescr',numpy.str_,256), #Description of measurements + ('Hdummy',numpy.str_,512), #Reserved space + #Main chunk + ('Msign','=5 + ('SPARrawGate2',' dateFile) or (endDate < dateFile): + continue + + self.fileList.append(thisFile) + self.dateFileList.append(dateFile) + + return + + def setNextFile(self): + + file_id = self.fileIndex + + if file_id == len(self.fileList): + print '\nNo more files in the folder' + print 'Total number of file(s) read : {}'.format(self.fileIndex + 1) + self.flagNoMoreFiles = 1 + return 0 + + print '\n[Setting file] (%s) ...' % self.fileList[file_id] + filename = os.path.join(self.path, self.fileList[file_id]) + + dirname, name = os.path.split(filename) + self.siteFile = name.split('.')[0] # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya + if self.filename is not None: + self.fp.close() + self.filename = filename + self.fp = open(self.filename, 'rb') + self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1) + self.nrecords = self.header_file['nrec'][0] + self.sizeOfFile = os.path.getsize(self.filename) + self.counter_records = 0 + self.flagIsNewFile = 0 + self.fileIndex += 1 + + return 1 + + def readNextBlock(self): + + while True: + if self.counter_records == self.nrecords: + self.flagIsNewFile = 1 + if not self.setNextFile(): + return 0 + + self.readBlock() + + if (self.datatime.time() < self.startTime) or (self.datatime.time() > self.endTime): + print "[Reading] Record No. %d/%d -> %s [Skipping]" %( + self.counter_records, + self.nrecords, + self.datatime.ctime()) + continue + break + + print "[Reading] Record No. %d/%d -> %s" %( + self.counter_records, + self.nrecords, + self.datatime.ctime()) + + return 1 + + def readBlock(self): + + pointer = self.fp.tell() + header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1) + self.nchannels = header_rec['nchan'][0]/2 + self.kchan = header_rec['nrxs'][0] + self.nmodes = header_rec['nmodes'][0] + self.nranges = header_rec['nranges'][0] + self.fp.seek(pointer) + self.height = numpy.empty((self.nmodes, self.nranges)) + self.snr = numpy.empty((self.nmodes, self.nchannels, self.nranges)) + self.buffer = numpy.empty((self.nmodes, 3, self.nranges)) + + for mode in range(self.nmodes): + self.readHeader() + data = self.readData() + self.height[mode] = (data[0] - self.correction) / 1000. + self.buffer[mode] = data[1] + self.snr[mode] = data[2] + + self.counter_records = self.counter_records + self.nmodes + + return + + def readHeader(self): + ''' + RecordHeader of BLTR rawdata file + ''' + + header_structure = numpy.dtype( + REC_HEADER_STRUCTURE.descr + [ + ('antenna_coord', 'f4', (2, self.nchannels)), + ('rx_gains', 'u4', (self.nchannels,)), + ('rx_analysis', 'u4', (self.nchannels,)) + ] + ) + + self.header_rec = numpy.fromfile(self.fp, header_structure, 1) + self.lat = self.header_rec['lat'][0] + self.lon = self.header_rec['lon'][0] + self.delta = self.header_rec['delta_r'][0] + self.correction = self.header_rec['dmode_rngcorr'][0] + self.imode = self.header_rec['dmode_index'][0] + self.antenna = self.header_rec['antenna_coord'] + self.rx_gains = self.header_rec['rx_gains'] + self.time = self.header_rec['time'][0] + tseconds = self.header_rec['time'][0] + local_t1 = time.localtime(tseconds) + self.year = local_t1.tm_year + self.month = local_t1.tm_mon + self.day = local_t1.tm_mday + self.t = datetime.datetime(self.year, self.month, self.day) + self.datatime = datetime.datetime.utcfromtimestamp(self.time) + + def readData(self): + ''' + Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value. + + Input: + status_value - Array data is set to NAN for values that are not equal to status_value + + ''' + + data_structure = numpy.dtype( + DATA_STRUCTURE.descr + [ + ('rx_saturation', 'u4', (self.nchannels,)), + ('chan_offset', 'u4', (2 * self.nchannels,)), + ('rx_amp', 'u4', (self.nchannels,)), + ('rx_snr', 'f4', (self.nchannels,)), + ('cross_snr', 'f4', (self.kchan,)), + ('sea_power_relative', 'f4', (self.kchan,))] + ) + + data = numpy.fromfile(self.fp, data_structure, self.nranges) + + height = data['range'] + winds = numpy.array((data['zonal'], data['meridional'], data['vertical'])) + snr = data['rx_snr'].T + + winds[numpy.where(winds == -9999.)] = numpy.nan + winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan + snr[numpy.where(snr == -9999.)] = numpy.nan + snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan + snr = numpy.power(10, snr / 10) + + return height, winds, snr + + def set_output(self): + ''' + Storing data from databuffer to dataOut object + ''' + + self.dataOut.data_SNR = self.snr + self.dataOut.height = self.height + self.dataOut.data_output = self.buffer + self.dataOut.utctimeInit = self.time + self.dataOut.utctime = self.dataOut.utctimeInit + self.dataOut.useLocalTime = False + self.dataOut.paramInterval = 157 + self.dataOut.timezone = self.timezone + self.dataOut.site = self.siteFile + self.dataOut.nrecords = self.nrecords/self.nmodes + self.dataOut.sizeOfFile = self.sizeOfFile + self.dataOut.lat = self.lat + self.dataOut.lon = self.lon + self.dataOut.channelList = range(self.nchannels) + self.dataOut.kchan = self.kchan + # self.dataOut.nHeights = self.nranges + self.dataOut.delta = self.delta + self.dataOut.correction = self.correction + self.dataOut.nmodes = self.nmodes + self.dataOut.imode = self.imode + self.dataOut.antenna = self.antenna + self.dataOut.rx_gains = self.rx_gains + self.dataOut.flagNoData = False + + def getData(self): + ''' + Storing data from databuffer to dataOut object + ''' + if self.flagNoMoreFiles: + self.dataOut.flagNoData = True + print 'No file left to process' + return 0 + + if not self.readNextBlock(): + self.dataOut.flagNoData = True + return 0 + + self.set_output() + + return 1 diff --git a/schainpy/model/io/jroIO_base.py b/schainpy/model/io/jroIO_base.py index d2c079c..ae5d5a1 100644 --- a/schainpy/model/io/jroIO_base.py +++ b/schainpy/model/io/jroIO_base.py @@ -1269,7 +1269,11 @@ class JRODataReader(JRODataIO): cursor=None, warnings=True, verbose=True, - server=None): + server=None, + format=None, + oneDDict=None, + twoDDict=None, + ind2DList=None): if server is not None: if 'tcp://' in server: address = server @@ -1417,11 +1421,12 @@ class JRODataReader(JRODataIO): print "[Reading] Number of read blocks %04d" %self.nTotalBlocks def printNumberOfBlock(self): + 'SPAM!' - if self.flagIsNewBlock: - print "[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks, - self.processingHeaderObj.dataBlocksPerFile, - self.dataOut.datatime.ctime()) +# if self.flagIsNewBlock: +# print "[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks, +# self.processingHeaderObj.dataBlocksPerFile, +# self.dataOut.datatime.ctime()) def printInfo(self): @@ -1456,7 +1461,11 @@ class JRODataReader(JRODataIO): cursor=None, warnings=True, server=None, - verbose=True, **kwargs): + verbose=True, + format=None, + oneDDict=None, + twoDDict=None, + ind2DList=None, **kwargs): if not(self.isConfig): self.setup(path=path, @@ -1479,7 +1488,11 @@ class JRODataReader(JRODataIO): cursor=cursor, warnings=warnings, server=server, - verbose=verbose) + verbose=verbose, + format=format, + oneDDict=oneDDict, + twoDDict=twoDDict, + ind2DList=ind2DList) self.isConfig = True if server is None: self.getData() diff --git a/schainpy/model/io/jroIO_bltr.py b/schainpy/model/io/jroIO_bltr.py new file mode 100644 index 0000000..fe46699 --- /dev/null +++ b/schainpy/model/io/jroIO_bltr.py @@ -0,0 +1,1154 @@ +import os, sys +import glob +import fnmatch +import datetime +import time +import re +import h5py +import numpy +import matplotlib.pyplot as plt + +import pylab as plb +from scipy.optimize import curve_fit +from scipy import asarray as ar, exp +from scipy import stats + +from numpy.ma.core import getdata + +SPEED_OF_LIGHT = 299792458 +SPEED_OF_LIGHT = 3e8 + +try: + from gevent import sleep +except: + from time import sleep + +from schainpy.model.data.jrodata import Spectra +#from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader +from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation +#from schainpy.model.io.jroIO_bltr import BLTRReader +from numpy import imag, shape, NaN + +from jroIO_base import JRODataReader + + +class Header(object): + + def __init__(self): + raise NotImplementedError + + + def read(self): + + raise NotImplementedError + + def write(self): + + raise NotImplementedError + + def printInfo(self): + + message = "#"*50 + "\n" + message += self.__class__.__name__.upper() + "\n" + message += "#"*50 + "\n" + + keyList = self.__dict__.keys() + keyList.sort() + + for key in keyList: + message += "%s = %s" %(key, self.__dict__[key]) + "\n" + + if "size" not in keyList: + attr = getattr(self, "size") + + if attr: + message += "%s = %s" %("size", attr) + "\n" + + #print message + + + + + +FILE_STRUCTURE = numpy.dtype([ #HEADER 48bytes + ('FileMgcNumber',' vertical) + ('AntennaCoord0',' endFp: + sys.stderr.write("Warning %s: Size value read from System Header is lower than it has to be\n" %fp) + return 0 + + if OffRHeader < endFp: + sys.stderr.write("Warning %s: Size value read from System Header size is greater than it has to be\n" %fp) + return 0 + + return 1 + + +class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODataReader): + + path = None + startDate = None + endDate = None + startTime = None + endTime = None + walk = None + isConfig = False + + + fileList= None + + #metadata + TimeZone= None + Interval= None + heightList= None + + #data + data= None + utctime= None + + + + def __init__(self, **kwargs): + + #Eliminar de la base la herencia + ProcessingUnit.__init__(self, **kwargs) + + #self.isConfig = False + + #self.pts2read_SelfSpectra = 0 + #self.pts2read_CrossSpectra = 0 + #self.pts2read_DCchannels = 0 + #self.datablock = None + self.utc = None + self.ext = ".fdt" + self.optchar = "P" + self.fpFile=None + self.fp = None + self.BlockCounter=0 + self.dtype = None + self.fileSizeByHeader = None + self.filenameList = [] + self.fileSelector = 0 + self.Off2StartNxtRec=0 + self.RecCounter=0 + self.flagNoMoreFiles = 0 + self.data_spc=None + self.data_cspc=None + self.data_output=None + self.path = None + self.OffsetStartHeader=0 + self.Off2StartData=0 + self.ipp = 0 + self.nFDTdataRecors=0 + self.blocksize = 0 + self.dataOut = Spectra() + self.profileIndex = 1 #Always + self.dataOut.flagNoData=False + self.dataOut.nRdPairs = 0 + self.dataOut.pairsList = [] + self.dataOut.data_spc=None + self.dataOut.noise=[] + self.dataOut.velocityX=[] + self.dataOut.velocityY=[] + self.dataOut.velocityV=[] + + + + def Files2Read(self, fp): + ''' + Function that indicates the number of .fdt files that exist in the folder to be read. + It also creates an organized list with the names of the files to read. + ''' + #self.__checkPath() + + ListaData=os.listdir(fp) #Gets the list of files within the fp address + ListaData=sorted(ListaData) #Sort the list of files from least to largest by names + nFiles=0 #File Counter + FileList=[] #A list is created that will contain the .fdt files + for IndexFile in ListaData : + if '.fdt' in IndexFile: + FileList.append(IndexFile) + nFiles+=1 + + #print 'Files2Read' + #print 'Existen '+str(nFiles)+' archivos .fdt' + + self.filenameList=FileList #List of files from least to largest by names + + + def run(self, **kwargs): + ''' + This method will be the one that will initiate the data entry, will be called constantly. + You should first verify that your Setup () is set up and then continue to acquire + the data to be processed with getData (). + ''' + if not self.isConfig: + self.setup(**kwargs) + self.isConfig = True + + self.getData() + #print 'running' + + + def setup(self, path=None, + startDate=None, + endDate=None, + startTime=None, + endTime=None, + walk=True, + timezone='utc', + code = None, + online=False, + ReadMode=None, + **kwargs): + + self.isConfig = True + + self.path=path + self.startDate=startDate + self.endDate=endDate + self.startTime=startTime + self.endTime=endTime + self.walk=walk + self.ReadMode=int(ReadMode) + + pass + + + def getData(self): + ''' + Before starting this function, you should check that there is still an unread file, + If there are still blocks to read or if the data block is empty. + + You should call the file "read". + + ''' + + if self.flagNoMoreFiles: + self.dataOut.flagNoData = True + print 'NoData se vuelve true' + return 0 + + self.fp=self.path + self.Files2Read(self.fp) + self.readFile(self.fp) + self.dataOut.data_spc = self.data_spc + self.dataOut.data_cspc =self.data_cspc + self.dataOut.data_output=self.data_output + + print 'self.dataOut.data_output', shape(self.dataOut.data_output) + + #self.removeDC() + return self.dataOut.data_spc + + + def readFile(self,fp): + ''' + You must indicate if you are reading in Online or Offline mode and load the + The parameters for this file reading mode. + + Then you must do 2 actions: + + 1. Get the BLTR FileHeader. + 2. Start reading the first block. + ''' + + #The address of the folder is generated the name of the .fdt file that will be read + print "File: ",self.fileSelector+1 + + if self.fileSelector < len(self.filenameList): + + self.fpFile=str(fp)+'/'+str(self.filenameList[self.fileSelector]) + #print self.fpFile + fheader = FileHeaderBLTR() + fheader.FHread(self.fpFile) #Bltr FileHeader Reading + self.nFDTdataRecors=fheader.nFDTdataRecors + + self.readBlock() #Block reading + else: + print 'readFile FlagNoData becomes true' + self.flagNoMoreFiles=True + self.dataOut.flagNoData = True + return 0 + + def getVelRange(self, extrapoints=0): + Lambda= SPEED_OF_LIGHT/50000000 + PRF = self.dataOut.PRF#1./(self.dataOut.ippSeconds * self.dataOut.nCohInt) + Vmax=-Lambda/(4.*(1./PRF)*self.dataOut.nCohInt*2.) + deltafreq = PRF / (self.nProfiles) + deltavel = (Vmax*2) / (self.nProfiles) + freqrange = deltafreq*(numpy.arange(self.nProfiles)-self.nProfiles/2.) - deltafreq/2 + velrange = deltavel*(numpy.arange(self.nProfiles)-self.nProfiles/2.) + return velrange + + def readBlock(self): + ''' + It should be checked if the block has data, if it is not passed to the next file. + + Then the following is done: + + 1. Read the RecordHeader + 2. Fill the buffer with the current block number. + + ''' + + if self.BlockCounter < self.nFDTdataRecors-2: + print self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + if self.ReadMode==1: + rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1) + elif self.ReadMode==0: + rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter) + + rheader.RHread(self.fpFile) #Bltr FileHeader Reading + + self.OffsetStartHeader=rheader.OffsetStartHeader + self.RecCounter=rheader.RecCounter + self.Off2StartNxtRec=rheader.Off2StartNxtRec + self.Off2StartData=rheader.Off2StartData + self.nProfiles=rheader.nProfiles + self.nChannels=rheader.nChannels + self.nHeights=rheader.nHeights + self.frequency=rheader.TransmitFrec + self.DualModeIndex=rheader.DualModeIndex + + self.pairsList =[(0,1),(0,2),(1,2)] + self.dataOut.pairsList = self.pairsList + + self.nRdPairs=len(self.dataOut.pairsList) + self.dataOut.nRdPairs = self.nRdPairs + + self.__firstHeigth=rheader.StartRangeSamp + self.__deltaHeigth=rheader.SampResolution + self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth + self.dataOut.channelList = range(self.nChannels) + self.dataOut.nProfiles=rheader.nProfiles + self.dataOut.nIncohInt=rheader.nIncohInt + self.dataOut.nCohInt=rheader.nCohInt + self.dataOut.ippSeconds= 1/float(rheader.PRFhz) + self.dataOut.PRF=rheader.PRFhz + self.dataOut.nFFTPoints=rheader.nProfiles + self.dataOut.utctime=rheader.nUtime + self.dataOut.timeZone=0 + self.dataOut.normFactor= self.dataOut.nProfiles*self.dataOut.nIncohInt*self.dataOut.nCohInt + self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles + + self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN + print 'self.data_output', shape(self.data_output) + self.dataOut.velocityX=[] + self.dataOut.velocityY=[] + self.dataOut.velocityV=[] + + '''Block Reading, the Block Data is received and Reshape is used to give it + shape. + ''' + + #Procedure to take the pointer to where the date block starts + startDATA = open(self.fpFile,"rb") + OffDATA= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec+self.Off2StartData + startDATA.seek(OffDATA, os.SEEK_SET) + + def moving_average(x, N=2): + return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):] + + def gaus(xSamples,a,x0,sigma): + return a*exp(-(xSamples-x0)**2/(2*sigma**2)) + + def Find(x,value): + for index in range(len(x)): + if x[index]==value: + return index + + def pol2cart(rho, phi): + x = rho * numpy.cos(phi) + y = rho * numpy.sin(phi) + return(x, y) + + + + + if self.DualModeIndex==self.ReadMode: + + self.data_fft = numpy.fromfile( startDATA, [('complex',' 0.0001) : +# +# try: +# popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) +# +# if numpy.amax(popt)>numpy.amax(yMean)*0.3: +# FitGauss=gaus(xSamples,*popt) +# +# else: +# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) +# print 'Verificador: Dentro', Height +# except RuntimeError: +# +# try: +# for j in range(len(ySamples[1])): +# yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]])) +# popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma]) +# FitGauss=gaus(xSamples,*popt) +# print 'Verificador: Exepcion1', Height +# except RuntimeError: +# +# try: +# popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma]) +# FitGauss=gaus(xSamples,*popt) +# print 'Verificador: Exepcion2', Height +# except RuntimeError: +# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) +# print 'Verificador: Exepcion3', Height +# else: +# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) +# #print 'Verificador: Fuera', Height +# +# +# +# Maximun=numpy.amax(yMean) +# eMinus1=Maximun*numpy.exp(-1) +# +# HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1))) +# HalfWidth= xFrec[HWpos] +# GCpos=Find(FitGauss, numpy.amax(FitGauss)) +# Vpos=Find(FactNorm, numpy.amax(FactNorm)) +# #Vpos=numpy.sum(FactNorm)/len(FactNorm) +# #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) ))) +# #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos +# '''****** Getting Fij ******''' +# +# GaussCenter=xFrec[GCpos] +# if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): +# Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 +# else: +# Fij=abs(GaussCenter-HalfWidth)+0.0000001 +# +# '''****** Getting Frecuency range of significant data ******''' +# +# Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10))) +# +# if Rangpos5 and len(FrecRange) 0.: +# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag +# #print 'Vmag',Vmag +# else: +# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN) +# +# if abs(Vx)<100 and abs(Vx) > 0.: +# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang +# #print 'Vang',Vang +# else: +# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN) +# +# if abs(GaussCenter)<2: +# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos]) +# +# else: +# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN) +# +# +# # print '********************************************' +# # print 'HalfWidth ', HalfWidth +# # print 'Maximun ', Maximun +# # print 'eMinus1 ', eMinus1 +# # print 'Rangpos ', Rangpos +# # print 'GaussCenter ',GaussCenter +# # print 'E01 ',E01 +# # print 'N01 ',N01 +# # print 'E02 ',E02 +# # print 'N02 ',N02 +# # print 'E12 ',E12 +# # print 'N12 ',N12 +# #print 'self.dataOut.velocityX ', self.dataOut.velocityX +# # print 'Fij ', Fij +# # print 'cC ', cC +# # print 'cF ', cF +# # print 'cG ', cG +# # print 'cA ', cA +# # print 'cB ', cB +# # print 'cH ', cH +# # print 'Vx ', Vx +# # print 'Vy ', Vy +# # print 'Vmag ', Vmag +# # print 'Vang ', Vang*180/numpy.pi +# # print 'PhaseSlope ',PhaseSlope[0] +# # print 'PhaseSlope ',PhaseSlope[1] +# # print 'PhaseSlope ',PhaseSlope[2] +# # print '********************************************' +# #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY) +# +# #print 'self.dataOut.velocityX', len(self.dataOut.velocityX) +# #print 'self.dataOut.velocityY', len(self.dataOut.velocityY) +# #print 'self.dataOut.velocityV', self.dataOut.velocityV +# +# self.data_output[0]=numpy.array(self.dataOut.velocityX) +# self.data_output[1]=numpy.array(self.dataOut.velocityY) +# self.data_output[2]=numpy.array(self.dataOut.velocityV) +# +# prin= self.data_output[0][~numpy.isnan(self.data_output[0])] +# print ' ' +# print 'VmagAverage',numpy.mean(prin) +# print ' ' +# # plt.figure(5) +# # plt.subplot(211) +# # plt.plot(self.dataOut.velocityX,'yo:') +# # plt.subplot(212) +# # plt.plot(self.dataOut.velocityY,'yo:') +# +# # plt.figure(1) +# # # plt.subplot(121) +# # # plt.plot(xFrec,ySamples[0],'k',label='Ch0') +# # # plt.plot(xFrec,ySamples[1],'g',label='Ch1') +# # # plt.plot(xFrec,ySamples[2],'r',label='Ch2') +# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') +# # # plt.legend() +# # plt.title('DATOS A ALTURA DE 2850 METROS') +# # +# # plt.xlabel('Frecuencia (KHz)') +# # plt.ylabel('Magnitud') +# # # plt.subplot(122) +# # # plt.title('Fit for Time Constant') +# # #plt.plot(xFrec,zline) +# # #plt.plot(xFrec,SmoothSPC,'g') +# # plt.plot(xFrec,FactNorm) +# # plt.axis([-4, 4, 0, 0.15]) +# # # plt.xlabel('SelfSpectra KHz') +# # +# # plt.figure(10) +# # # plt.subplot(121) +# # plt.plot(xFrec,ySamples[0],'b',label='Ch0') +# # plt.plot(xFrec,ySamples[1],'y',label='Ch1') +# # plt.plot(xFrec,ySamples[2],'r',label='Ch2') +# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') +# # plt.legend() +# # plt.title('SELFSPECTRA EN CANALES') +# # +# # plt.xlabel('Frecuencia (KHz)') +# # plt.ylabel('Magnitud') +# # # plt.subplot(122) +# # # plt.title('Fit for Time Constant') +# # #plt.plot(xFrec,zline) +# # #plt.plot(xFrec,SmoothSPC,'g') +# # # plt.plot(xFrec,FactNorm) +# # # plt.axis([-4, 4, 0, 0.15]) +# # # plt.xlabel('SelfSpectra KHz') +# # +# # plt.figure(9) +# # +# # +# # plt.title('DATOS SUAVIZADOS') +# # plt.xlabel('Frecuencia (KHz)') +# # plt.ylabel('Magnitud') +# # plt.plot(xFrec,SmoothSPC,'g') +# # +# # #plt.plot(xFrec,FactNorm) +# # plt.axis([-4, 4, 0, 0.15]) +# # # plt.xlabel('SelfSpectra KHz') +# # # +# # plt.figure(2) +# # # #plt.subplot(121) +# # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra') +# # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano') +# # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo') +# # # #plt.plot(xFrec,phase) +# # # plt.xlabel('Suavizado, promediado KHz') +# # plt.title('SELFSPECTRA PROMEDIADO') +# # # #plt.subplot(122) +# # # #plt.plot(xSamples,zline) +# # plt.xlabel('Frecuencia (KHz)') +# # plt.ylabel('Magnitud') +# # plt.legend() +# # # +# # # plt.figure(3) +# # # plt.subplot(311) +# # # #plt.plot(xFrec,phase[0]) +# # # plt.plot(xFrec,phase[0],'g') +# # # plt.subplot(312) +# # # plt.plot(xFrec,phase[1],'g') +# # # plt.subplot(313) +# # # plt.plot(xFrec,phase[2],'g') +# # # #plt.plot(xFrec,phase[2]) +# # # +# # # plt.figure(4) +# # # +# # # plt.plot(xSamples,coherence[0],'b') +# # # plt.plot(xSamples,coherence[1],'r') +# # # plt.plot(xSamples,coherence[2],'g') +# # plt.show() +# # # +# # # plt.clf() +# # # plt.cla() +# # # plt.close() +# +# print ' ' + + + + self.BlockCounter+=2 + + else: + self.fileSelector+=1 + self.BlockCounter=0 + print "Next File" + + + + + diff --git a/schainpy/model/io/jroIO_madrigal.py b/schainpy/model/io/jroIO_madrigal.py new file mode 100644 index 0000000..7a62ffe --- /dev/null +++ b/schainpy/model/io/jroIO_madrigal.py @@ -0,0 +1,632 @@ +''' +Created on Aug 1, 2017 + +@author: Juan C. Espinoza +''' + +import os +import sys +import time +import json +import glob +import datetime + +import numpy +import h5py + +from schainpy.model.io.jroIO_base import JRODataReader +from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation +from schainpy.model.data.jrodata import Parameters +from schainpy.utils import log + +try: + import madrigal.cedar +except: + log.warning( + 'You should install "madrigal library" module if you want to read/write Madrigal data' + ) + +DEF_CATALOG = { + 'principleInvestigator': 'Marco Milla', + 'expPurpose': None, + 'cycleTime': None, + 'correlativeExp': None, + 'sciRemarks': None, + 'instRemarks': None + } +DEF_HEADER = { + 'kindatDesc': None, + 'analyst': 'Jicamarca User', + 'comments': None, + 'history': None + } +MNEMONICS = { + 10: 'jro', + 11: 'jbr', + 840: 'jul', + 13: 'jas', + 1000: 'pbr', + 1001: 'hbr', + 1002: 'obr', +} + +UT1970 = datetime.datetime(1970, 1, 1) - datetime.timedelta(seconds=time.timezone) + +def load_json(obj): + ''' + Parse json as string instead of unicode + ''' + + if isinstance(obj, str): + iterable = json.loads(obj) + else: + iterable = obj + + if isinstance(iterable, dict): + return {str(k): load_json(v) if isinstance(v, dict) else str(v) if isinstance(v, unicode) else v + for k, v in iterable.items()} + elif isinstance(iterable, (list, tuple)): + return [str(v) if isinstance(v, unicode) else v for v in iterable] + + return iterable + + +class MADReader(JRODataReader, ProcessingUnit): + + def __init__(self, **kwargs): + + ProcessingUnit.__init__(self, **kwargs) + + self.dataOut = Parameters() + self.counter_records = 0 + self.nrecords = None + self.flagNoMoreFiles = 0 + self.isConfig = False + self.filename = None + self.intervals = set() + + def setup(self, + path=None, + startDate=None, + endDate=None, + format=None, + startTime=datetime.time(0, 0, 0), + endTime=datetime.time(23, 59, 59), + **kwargs): + + self.path = path + self.startDate = startDate + self.endDate = endDate + self.startTime = startTime + self.endTime = endTime + self.datatime = datetime.datetime(1900,1,1) + self.oneDDict = load_json(kwargs.get('oneDDict', + "{\"GDLATR\":\"lat\", \"GDLONR\":\"lon\"}")) + self.twoDDict = load_json(kwargs.get('twoDDict', + "{\"GDALT\": \"heightList\"}")) + self.ind2DList = load_json(kwargs.get('ind2DList', + "[\"GDALT\"]")) + if self.path is None: + raise ValueError, 'The path is not valid' + + if format is None: + raise ValueError, 'The format is not valid choose simple or hdf5' + elif format.lower() in ('simple', 'txt'): + self.ext = '.txt' + elif format.lower() in ('cedar',): + self.ext = '.001' + else: + self.ext = '.hdf5' + + self.search_files(self.path) + self.fileId = 0 + + if not self.fileList: + raise Warning, 'There is no files matching these date in the folder: {}. \n Check startDate and endDate'.format(path) + + self.setNextFile() + + def search_files(self, path): + ''' + Searching for madrigal files in path + Creating a list of files to procces included in [startDate,endDate] + + Input: + path - Path to find files + ''' + + log.log('Searching files {} in {} '.format(self.ext, path), 'MADReader') + foldercounter = 0 + fileList0 = glob.glob1(path, '*{}'.format(self.ext)) + fileList0.sort() + + self.fileList = [] + self.dateFileList = [] + + startDate = self.startDate - datetime.timedelta(1) + endDate = self.endDate + datetime.timedelta(1) + + for thisFile in fileList0: + year = thisFile[3:7] + if not year.isdigit(): + continue + + month = thisFile[7:9] + if not month.isdigit(): + continue + + day = thisFile[9:11] + if not day.isdigit(): + continue + + year, month, day = int(year), int(month), int(day) + dateFile = datetime.date(year, month, day) + + if (startDate > dateFile) or (endDate < dateFile): + continue + + self.fileList.append(thisFile) + self.dateFileList.append(dateFile) + + return + + def parseHeader(self): + ''' + ''' + + self.output = {} + self.version = '2' + s_parameters = None + if self.ext == '.txt': + self.parameters = [s.strip().lower() for s in self.fp.readline().strip().split(' ') if s] + elif self.ext == '.hdf5': + metadata = self.fp['Metadata'] + data = self.fp['Data']['Array Layout'] + if 'Independent Spatial Parameters' in metadata: + s_parameters = [s[0].lower() for s in metadata['Independent Spatial Parameters']] + self.version = '3' + one = [s[0].lower() for s in data['1D Parameters']['Data Parameters']] + one_d = [1 for s in one] + two = [s[0].lower() for s in data['2D Parameters']['Data Parameters']] + two_d = [2 for s in two] + self.parameters = one + two + self.parameters_d = one_d + two_d + + log.success('Parameters found: {}'.format(','.join(self.parameters)), + 'MADReader') + if s_parameters: + log.success('Spatial parameters: {}'.format(','.join(s_parameters)), + 'MADReader') + + for param in self.oneDDict.keys(): + if param.lower() not in self.parameters: + log.warning( + 'Parameter {} not found will be ignored'.format( + param), + 'MADReader') + self.oneDDict.pop(param, None) + + for param, value in self.twoDDict.items(): + if param.lower() not in self.parameters: + log.warning( + 'Parameter {} not found, it will be ignored'.format( + param), + 'MADReader') + self.twoDDict.pop(param, None) + continue + if isinstance(value, list): + if value[0] not in self.output: + self.output[value[0]] = [] + self.output[value[0]].append(None) + + def parseData(self): + ''' + ''' + + if self.ext == '.txt': + self.data = numpy.genfromtxt(self.fp, missing_values=('missing')) + self.nrecords = self.data.shape[0] + self.ranges = numpy.unique(self.data[:,self.parameters.index(self.ind2DList[0].lower())]) + elif self.ext == '.hdf5': + self.data = self.fp['Data']['Array Layout'] + self.nrecords = len(self.data['timestamps'].value) + self.ranges = self.data['range'].value + + def setNextFile(self): + ''' + ''' + + file_id = self.fileId + + if file_id == len(self.fileList): + log.success('No more files', 'MADReader') + self.flagNoMoreFiles = 1 + return 0 + + log.success( + 'Opening: {}'.format(self.fileList[file_id]), + 'MADReader' + ) + + filename = os.path.join(self.path, self.fileList[file_id]) + + if self.filename is not None: + self.fp.close() + + self.filename = filename + self.filedate = self.dateFileList[file_id] + + if self.ext=='.hdf5': + self.fp = h5py.File(self.filename, 'r') + else: + self.fp = open(self.filename, 'rb') + + self.parseHeader() + self.parseData() + self.sizeOfFile = os.path.getsize(self.filename) + self.counter_records = 0 + self.flagIsNewFile = 0 + self.fileId += 1 + + return 1 + + def readNextBlock(self): + + while True: + self.flagDiscontinuousBlock = 0 + if self.flagIsNewFile: + if not self.setNextFile(): + return 0 + + self.readBlock() + + if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \ + (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)): + log.warning( + 'Reading Record No. {}/{} -> {} [Skipping]'.format( + self.counter_records, + self.nrecords, + self.datatime.ctime()), + 'MADReader') + continue + break + + log.log( + 'Reading Record No. {}/{} -> {}'.format( + self.counter_records, + self.nrecords, + self.datatime.ctime()), + 'MADReader') + + return 1 + + def readBlock(self): + ''' + ''' + dum = [] + if self.ext == '.txt': + dt = self.data[self.counter_records][:6].astype(int) + self.datatime = datetime.datetime(dt[0], dt[1], dt[2], dt[3], dt[4], dt[5]) + while True: + dt = self.data[self.counter_records][:6].astype(int) + datatime = datetime.datetime(dt[0], dt[1], dt[2], dt[3], dt[4], dt[5]) + if datatime == self.datatime: + dum.append(self.data[self.counter_records]) + self.counter_records += 1 + if self.counter_records == self.nrecords: + self.flagIsNewFile = True + break + continue + self.intervals.add((datatime-self.datatime).seconds) + if datatime.date() > self.datatime.date(): + self.flagDiscontinuousBlock = 1 + break + elif self.ext == '.hdf5': + datatime = datetime.datetime.utcfromtimestamp( + self.data['timestamps'][self.counter_records]) + nHeights = len(self.ranges) + for n, param in enumerate(self.parameters): + if self.parameters_d[n] == 1: + dum.append(numpy.ones(nHeights)*self.data['1D Parameters'][param][self.counter_records]) + else: + if self.version == '2': + dum.append(self.data['2D Parameters'][param][self.counter_records]) + else: + tmp = self.data['2D Parameters'][param].value.T + dum.append(tmp[self.counter_records]) + self.intervals.add((datatime-self.datatime).seconds) + if datatime.date()>self.datatime.date(): + self.flagDiscontinuousBlock = 1 + self.datatime = datatime + self.counter_records += 1 + if self.counter_records == self.nrecords: + self.flagIsNewFile = True + + self.buffer = numpy.array(dum) + return + + def set_output(self): + ''' + Storing data from buffer to dataOut object + ''' + + parameters = [None for __ in self.parameters] + + for param, attr in self.oneDDict.items(): + x = self.parameters.index(param.lower()) + setattr(self.dataOut, attr, self.buffer[0][x]) + + for param, value in self.twoDDict.items(): + x = self.parameters.index(param.lower()) + if self.ext == '.txt': + y = self.parameters.index(self.ind2DList[0].lower()) + ranges = self.buffer[:,y] + if self.ranges.size == ranges.size: + continue + index = numpy.where(numpy.in1d(self.ranges, ranges))[0] + dummy = numpy.zeros(self.ranges.shape) + numpy.nan + dummy[index] = self.buffer[:,x] + else: + dummy = self.buffer[x] + + if isinstance(value, str): + if value not in self.ind2DList: + setattr(self.dataOut, value, dummy.reshape(1,-1)) + elif isinstance(value, list): + self.output[value[0]][value[1]] = dummy + parameters[value[1]] = param + + for key, value in self.output.items(): + setattr(self.dataOut, key, numpy.array(value)) + + self.dataOut.parameters = [s for s in parameters if s] + self.dataOut.heightList = self.ranges + self.dataOut.utctime = (self.datatime - UT1970).total_seconds() + self.dataOut.utctimeInit = self.dataOut.utctime + self.dataOut.paramInterval = min(self.intervals) + self.dataOut.useLocalTime = False + self.dataOut.flagNoData = False + self.dataOut.nrecords = self.nrecords + self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock + + def getData(self): + ''' + Storing data from databuffer to dataOut object + ''' + if self.flagNoMoreFiles: + self.dataOut.flagNoData = True + log.error('No file left to process', 'MADReader') + return 0 + + if not self.readNextBlock(): + self.dataOut.flagNoData = True + return 0 + + self.set_output() + + return 1 + + +class MADWriter(Operation): + + missing = -32767 + + def __init__(self, **kwargs): + + Operation.__init__(self, **kwargs) + self.dataOut = Parameters() + self.path = None + self.fp = None + + def run(self, dataOut, path, oneDDict, ind2DList='[]', twoDDict='{}', + metadata='{}', format='cedar', **kwargs): + ''' + Inputs: + path - path where files will be created + oneDDict - json of one-dimensional parameters in record where keys + are Madrigal codes (integers or mnemonics) and values the corresponding + dataOut attribute e.g: { + 'gdlatr': 'lat', + 'gdlonr': 'lon', + 'gdlat2':'lat', + 'glon2':'lon'} + ind2DList - list of independent spatial two-dimensional parameters e.g: + ['heighList'] + twoDDict - json of two-dimensional parameters in record where keys + are Madrigal codes (integers or mnemonics) and values the corresponding + dataOut attribute if multidimensional array specify as tupple + ('attr', pos) e.g: { + 'gdalt': 'heightList', + 'vn1p2': ('data_output', 0), + 'vn2p2': ('data_output', 1), + 'vn3': ('data_output', 2), + 'snl': ('data_SNR', 'db') + } + metadata - json of madrigal metadata (kinst, kindat, catalog and header) + ''' + if not self.isConfig: + self.setup(path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs) + self.isConfig = True + + self.dataOut = dataOut + self.putData() + return + + def setup(self, path, oneDDict, ind2DList, twoDDict, metadata, format, **kwargs): + ''' + Configure Operation + ''' + + self.path = path + self.blocks = kwargs.get('blocks', None) + self.counter = 0 + self.oneDDict = load_json(oneDDict) + self.twoDDict = load_json(twoDDict) + self.ind2DList = load_json(ind2DList) + meta = load_json(metadata) + self.kinst = meta.get('kinst') + self.kindat = meta.get('kindat') + self.catalog = meta.get('catalog', DEF_CATALOG) + self.header = meta.get('header', DEF_HEADER) + if format == 'cedar': + self.ext = '.dat' + self.extra_args = {} + elif format == 'hdf5': + self.ext = '.hdf5' + self.extra_args = {'ind2DList': self.ind2DList} + + self.keys = [k.lower() for k in self.twoDDict] + if 'range' in self.keys: + self.keys.remove('range') + if 'gdalt' in self.keys: + self.keys.remove('gdalt') + + def setFile(self): + ''' + Create new cedar file object + ''' + + self.mnemonic = MNEMONICS[self.kinst] #TODO get mnemonic from madrigal + date = datetime.datetime.fromtimestamp(self.dataOut.utctime) + + filename = '{}{}{}'.format(self.mnemonic, + date.strftime('%Y%m%d_%H%M%S'), + self.ext) + + self.fullname = os.path.join(self.path, filename) + + if os.path.isfile(self.fullname) : + log.warning( + 'Destination path {} already exists. Previous file deleted.'.format( + self.fullname), + 'MADWriter') + os.remove(self.fullname) + + try: + log.success( + 'Creating file: {}'.format(self.fullname), + 'MADWriter') + self.fp = madrigal.cedar.MadrigalCedarFile(self.fullname, True) + except ValueError, e: + log.error( + 'Impossible to create a cedar object with "madrigal.cedar.MadrigalCedarFile"', + 'MADWriter') + return + + return 1 + + def writeBlock(self): + ''' + Add data records to cedar file taking data from oneDDict and twoDDict + attributes. + Allowed parameters in: parcodes.tab + ''' + + startTime = datetime.datetime.fromtimestamp(self.dataOut.utctime) + endTime = startTime + datetime.timedelta(seconds=self.dataOut.paramInterval) + heights = self.dataOut.heightList + + if self.ext == '.dat': + invalid = numpy.isnan(self.dataOut.data_output) + self.dataOut.data_output[invalid] = self.missing + out = {} + for key, value in self.twoDDict.items(): + key = key.lower() + if isinstance(value, str): + if 'db' in value.lower(): + tmp = getattr(self.dataOut, value.replace('_db', '')) + SNRavg = numpy.average(tmp, axis=0) + tmp = 10*numpy.log10(SNRavg) + else: + tmp = getattr(self.dataOut, value) + out[key] = tmp.flatten() + elif isinstance(value, (tuple, list)): + attr, x = value + data = getattr(self.dataOut, attr) + out[key] = data[int(x)] + + a = numpy.array([out[k] for k in self.keys]) + nrows = numpy.array([numpy.isnan(a[:, x]).all() for x in range(len(heights))]) + index = numpy.where(nrows == False)[0] + + rec = madrigal.cedar.MadrigalDataRecord( + self.kinst, + self.kindat, + startTime.year, + startTime.month, + startTime.day, + startTime.hour, + startTime.minute, + startTime.second, + startTime.microsecond/10000, + endTime.year, + endTime.month, + endTime.day, + endTime.hour, + endTime.minute, + endTime.second, + endTime.microsecond/10000, + self.oneDDict.keys(), + self.twoDDict.keys(), + len(index), + **self.extra_args + ) + + # Setting 1d values + for key in self.oneDDict: + rec.set1D(key, getattr(self.dataOut, self.oneDDict[key])) + + # Setting 2d values + nrec = 0 + for n in index: + for key in out: + rec.set2D(key, nrec, out[key][n]) + nrec += 1 + + self.fp.append(rec) + if self.ext == '.hdf5' and self.counter % 500 == 0 and self.counter > 0: + self.fp.dump() + if self.counter % 10 == 0 and self.counter > 0: + log.log( + 'Writing {} records'.format( + self.counter), + 'MADWriter') + + def setHeader(self): + ''' + Create an add catalog and header to cedar file + ''' + + log.success('Closing file {}'.format(self.fullname), 'MADWriter') + + if self.ext == '.dat': + self.fp.write() + else: + self.fp.dump() + self.fp.close() + + header = madrigal.cedar.CatalogHeaderCreator(self.fullname) + header.createCatalog(**self.catalog) + header.createHeader(**self.header) + header.write() + + def putData(self): + + if self.dataOut.flagNoData: + return 0 + + if self.dataOut.flagDiscontinuousBlock or self.counter == self.blocks: + if self.counter > 0: + self.setHeader() + self.counter = 0 + + if self.counter == 0: + self.setFile() + + self.writeBlock() + self.counter += 1 + + def close(self): + + if self.counter > 0: + self.setHeader() diff --git a/schainpy/model/io/jroIO_mira35c.py b/schainpy/model/io/jroIO_mira35c.py new file mode 100644 index 0000000..632bc93 --- /dev/null +++ b/schainpy/model/io/jroIO_mira35c.py @@ -0,0 +1,803 @@ +import os, sys +import glob +import fnmatch +import datetime +import time +import re +import h5py +import numpy +import matplotlib.pyplot as plt + +import pylab as plb +from scipy.optimize import curve_fit +from scipy import asarray as ar,exp +from scipy import stats + +from numpy.ma.core import getdata + +SPEED_OF_LIGHT = 299792458 +SPEED_OF_LIGHT = 3e8 + +try: + from gevent import sleep +except: + from time import sleep + +from schainpy.model.data.jrodata import Spectra +#from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader +from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation +#from schainpy.model.io.jroIO_bltr import BLTRReader +from numpy import imag, shape, NaN, empty + + + +class Header(object): + + def __init__(self): + raise NotImplementedError + + + def read(self): + + raise NotImplementedError + + def write(self): + + raise NotImplementedError + + def printInfo(self): + + message = "#"*50 + "\n" + message += self.__class__.__name__.upper() + "\n" + message += "#"*50 + "\n" + + keyList = self.__dict__.keys() + keyList.sort() + + for key in keyList: + message += "%s = %s" %(key, self.__dict__[key]) + "\n" + + if "size" not in keyList: + attr = getattr(self, "size") + + if attr: + message += "%s = %s" %("size", attr) + "\n" + + #print message + + +FILE_HEADER = numpy.dtype([ #HEADER 1024bytes + ('Hname','a32'), #Original file name + ('Htime',numpy.str_,32), #Date and time when the file was created + ('Hoper',numpy.str_,64), #Name of operator who created the file + ('Hplace',numpy.str_,128), #Place where the measurements was carried out + ('Hdescr',numpy.str_,256), #Description of measurements + ('Hdummy',numpy.str_,512), #Reserved space + #Main chunk 8bytes + ('Msign',numpy.str_,4), #Main chunk signature FZKF or NUIG + ('MsizeData','=5 + ('SPARrawGate2',' 1180: + self.fp.seek(self.PointerReader , os.SEEK_SET) + self.FirstPoint = self.PointerReader + + else : + self.FirstPoint = 1180 + + + + self.srviHeader = SRVIHeader() + + self.srviHeader.SRVIread(self.fp) #Se obtiene la cabecera del SRVI + + self.blocksize = self.srviHeader.SizeOfDataBlock1 # Se obtiene el tamao del bloque + + if self.blocksize == 148: + print 'blocksize == 148 bug' + jump = numpy.fromfile(self.fp,[('jump',numpy.str_,140)] ,1) + + self.srviHeader.SRVIread(self.fp) #Se obtiene la cabecera del SRVI + + if not self.srviHeader.SizeOfSRVI1: + self.fileSelector+=1 + self.nextfileflag==True + self.FileHeaderFlag == True + + self.recordheader = RecordHeader() + self.recordheader.RHread(self.fp) + self.RadarConst = self.recordheader.RadarConst + dwell = self.recordheader.time_t + npw1 = self.recordheader.npw1 + npw2 = self.recordheader.npw2 + + + self.dataOut.channelList = range(1) + self.dataOut.nIncohInt = self.Num_inCoh + self.dataOut.nProfiles = self.Num_Bins + self.dataOut.nCohInt = 1 + self.dataOut.windowOfFilter = 1 + self.dataOut.utctime = dwell + self.dataOut.timeZone=0 + + self.dataOut.outputInterval = self.dataOut.getTimeInterval() + self.dataOut.heightList = self.SPARrawGate1*self.__deltaHeigth + numpy.array(range(self.Num_Hei))*self.__deltaHeigth + + + + self.HSDVsign = numpy.fromfile( self.fp, [('HSDV',numpy.str_,4)],1) + self.SizeHSDV = numpy.fromfile( self.fp, [('SizeHSDV',' 0. , self.data_spc, 0) + + self.dataOut.COFA = numpy.array([self.COFA_Co , self.COFA_Cx]) + + print ' ' + print 'SPC',numpy.shape(self.dataOut.data_spc) + #print 'SPC',self.dataOut.data_spc + + noinor1 = 713031680 + noinor2 = 30 + + npw1 = 1#0**(npw1/10) * noinor1 * noinor2 + npw2 = 1#0**(npw2/10) * noinor1 * noinor2 + self.dataOut.NPW = numpy.array([npw1, npw2]) + + print ' ' + + self.data_spc = numpy.transpose(self.data_spc, (2,1,0)) + self.data_spc = numpy.fft.fftshift(self.data_spc, axes = 1) + + self.data_spc = numpy.fliplr(self.data_spc) + + self.data_spc = numpy.where(self.data_spc > 0. , self.data_spc, 0) + self.dataOut_spc= numpy.ones([1, self.Num_Bins , self.Num_Hei]) + self.dataOut_spc[0,:,:] = self.data_spc[0,:,:] + #print 'SHAPE', self.dataOut_spc.shape + #For nyquist correction: + #fix = 20 # ~3m/s + #shift = self.Num_Bins/2 + fix + #self.data_spc = numpy.array([ self.data_spc[: , self.Num_Bins-shift+1: , :] , self.data_spc[: , 0:self.Num_Bins-shift , :]]) + + + + '''Block Reading, the Block Data is received and Reshape is used to give it + shape. + ''' + + self.PointerReader = self.fp.tell() + + + + + + + \ No newline at end of file diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index 50e7367..fa4495c 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -178,8 +178,8 @@ class ParamReader(ProcessingUnit): print "[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime) print - for i in range(len(filenameList)): - print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime()) +# for i in range(len(filenameList)): +# print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime()) self.filenameList = filenameList self.datetimeList = datetimeList diff --git a/schainpy/model/io/jroIO_spectra.py b/schainpy/model/io/jroIO_spectra.py index d4f55df..9daf643 100644 --- a/schainpy/model/io/jroIO_spectra.py +++ b/schainpy/model/io/jroIO_spectra.py @@ -11,6 +11,7 @@ from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, from schainpy.model.data.jrodata import Spectra class SpectraReader(JRODataReader, ProcessingUnit): + """ Esta clase permite leer datos de espectros desde archivos procesados (.pdata). La lectura de los datos siempre se realiza por bloques. Los datos leidos (array de 3 dimensiones) @@ -20,6 +21,7 @@ class SpectraReader(JRODataReader, ProcessingUnit): paresCanalesDiferentes * alturas * perfiles (Cross Spectra) canales * alturas (DC Channels) + Esta clase contiene instancias (objetos) de las clases BasicHeader, SystemHeader, RadarControllerHeader y Spectra. Los tres primeros se usan para almacenar informacion de la cabecera de datos (metadata), y el cuarto (Spectra) para obtener y almacenar un bloque de @@ -74,6 +76,7 @@ class SpectraReader(JRODataReader, ProcessingUnit): Inicializador de la clase SpectraReader para la lectura de datos de espectros. Inputs: + dataOut : Objeto de la clase Spectra. Este objeto sera utilizado para almacenar un perfil de datos cada vez que se haga un requerimiento (getData). El perfil sera obtenido a partir del buffer de datos, @@ -81,104 +84,107 @@ class SpectraReader(JRODataReader, ProcessingUnit): bloque de datos. Si este parametro no es pasado se creara uno internamente. - Affected: + + Affected: + self.dataOut Return : None """ + #Eliminar de la base la herencia ProcessingUnit.__init__(self, **kwargs) - + # self.isConfig = False - + self.pts2read_SelfSpectra = 0 - + self.pts2read_CrossSpectra = 0 - + self.pts2read_DCchannels = 0 - + self.datablock = None - + self.utc = None - + self.ext = ".pdata" - + self.optchar = "P" - + self.basicHeaderObj = BasicHeader(LOCALTIME) - + self.systemHeaderObj = SystemHeader() - + self.radarControllerHeaderObj = RadarControllerHeader() - + self.processingHeaderObj = ProcessingHeader() - + self.online = 0 - + self.fp = None - + self.idFile = None - + self.dtype = None - + self.fileSizeByHeader = None - + self.filenameList = [] - + self.filename = None - + self.fileSize = None - + self.firstHeaderSize = 0 - + self.basicHeaderSize = 24 - + self.pathList = [] self.lastUTTime = 0 - + self.maxTimeStep = 30 - + self.flagNoMoreFiles = 0 - + self.set = 0 - + self.path = None self.delay = 60 #seconds - + self.nTries = 3 #quantity tries - + self.nFiles = 3 #number of files for searching - + self.nReadBlocks = 0 - + self.flagIsNewFile = 1 - + self.__isFirstTimeOnline = 1 - + # self.ippSeconds = 0 - - self.flagDiscontinuousBlock = 0 - + + self.flagDiscontinuousBlock = 0 + self.flagIsNewBlock = 0 - + self.nTotalBlocks = 0 - + self.blocksize = 0 - + self.dataOut = self.createObjByDefault() - + self.profileIndex = 1 #Always def createObjByDefault(self): - + dataObj = Spectra() - + return dataObj - + def __hasNotDataInBuffer(self): return 1 @@ -186,7 +192,7 @@ class SpectraReader(JRODataReader, ProcessingUnit): def getBlockDimension(self): """ Obtiene la cantidad de puntos a leer por cada bloque de datos - + Affected: self.nRdChannels self.nRdPairs @@ -206,37 +212,39 @@ class SpectraReader(JRODataReader, ProcessingUnit): for i in range(0, self.processingHeaderObj.totalSpectra*2, 2): if self.processingHeaderObj.spectraComb[i] == self.processingHeaderObj.spectraComb[i+1]: - self.nRdChannels = self.nRdChannels + 1 #par de canales iguales + self.nRdChannels = self.nRdChannels + 1 #par de canales iguales + else: self.nRdPairs = self.nRdPairs + 1 #par de canales diferentes self.rdPairList.append((self.processingHeaderObj.spectraComb[i], self.processingHeaderObj.spectraComb[i+1])) pts2read = self.processingHeaderObj.nHeights * self.processingHeaderObj.profilesPerBlock - + self.pts2read_SelfSpectra = int(self.nRdChannels * pts2read) self.blocksize = self.pts2read_SelfSpectra - + if self.processingHeaderObj.flag_cspc: self.pts2read_CrossSpectra = int(self.nRdPairs * pts2read) self.blocksize += self.pts2read_CrossSpectra - + if self.processingHeaderObj.flag_dc: self.pts2read_DCchannels = int(self.systemHeaderObj.nChannels * self.processingHeaderObj.nHeights) self.blocksize += self.pts2read_DCchannels - + # self.blocksize = self.pts2read_SelfSpectra + self.pts2read_CrossSpectra + self.pts2read_DCchannels - + def readBlock(self): """ Lee el bloque de datos desde la posicion actual del puntero del archivo (self.fp) y actualiza todos los parametros relacionados al bloque de datos (metadata + data). La data leida es almacenada en el buffer y el contador del buffer es seteado a 0 - + Return: None - + Variables afectadas: + self.flagIsNewFile self.flagIsNewBlock @@ -245,9 +253,10 @@ class SpectraReader(JRODataReader, ProcessingUnit): self.data_cspc self.data_dc - Exceptions: + Exceptions: Si un bloque leido no es un bloque valido """ + blockOk_flag = False fpointer = self.fp.tell() @@ -257,30 +266,32 @@ class SpectraReader(JRODataReader, ProcessingUnit): if self.processingHeaderObj.flag_cspc: cspc = numpy.fromfile( self.fp, self.dtype, self.pts2read_CrossSpectra ) cspc = cspc.reshape( (self.nRdPairs, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D - + if self.processingHeaderObj.flag_dc: dc = numpy.fromfile( self.fp, self.dtype, self.pts2read_DCchannels ) #int(self.processingHeaderObj.nHeights*self.systemHeaderObj.nChannels) ) dc = dc.reshape( (self.systemHeaderObj.nChannels, self.processingHeaderObj.nHeights) ) #transforma a un arreglo 2D - - + + if not(self.processingHeaderObj.shif_fft): #desplaza a la derecha en el eje 2 determinadas posiciones shift = int(self.processingHeaderObj.profilesPerBlock/2) spc = numpy.roll( spc, shift , axis=2 ) - + if self.processingHeaderObj.flag_cspc: #desplaza a la derecha en el eje 2 determinadas posiciones cspc = numpy.roll( cspc, shift, axis=2 ) - + #Dimensions : nChannels, nProfiles, nSamples spc = numpy.transpose( spc, (0,2,1) ) self.data_spc = spc + + if self.processingHeaderObj.flag_cspc: - if self.processingHeaderObj.flag_cspc: cspc = numpy.transpose( cspc, (0,2,1) ) self.data_cspc = cspc['real'] + cspc['imag']*1j else: self.data_cspc = None + if self.processingHeaderObj.flag_dc: self.data_dc = dc['real'] + dc['imag']*1j @@ -294,60 +305,60 @@ class SpectraReader(JRODataReader, ProcessingUnit): self.nReadBlocks += 1 return 1 - + def getFirstHeader(self): - + self.getBasicHeader() - + self.dataOut.systemHeaderObj = self.systemHeaderObj.copy() - + self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy() - + # self.dataOut.ippSeconds = self.ippSeconds - + # self.dataOut.timeInterval = self.radarControllerHeaderObj.ippSeconds * self.processingHeaderObj.nCohInt * self.processingHeaderObj.nIncohInt * self.processingHeaderObj.profilesPerBlock self.dataOut.dtype = self.dtype - + # self.dataOut.nPairs = self.nPairs - + self.dataOut.pairsList = self.rdPairList - + self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock - + self.dataOut.nFFTPoints = self.processingHeaderObj.profilesPerBlock - + self.dataOut.nCohInt = self.processingHeaderObj.nCohInt - + self.dataOut.nIncohInt = self.processingHeaderObj.nIncohInt - + xf = self.processingHeaderObj.firstHeight + self.processingHeaderObj.nHeights*self.processingHeaderObj.deltaHeight - self.dataOut.heightList = numpy.arange(self.processingHeaderObj.firstHeight, xf, self.processingHeaderObj.deltaHeight) - + self.dataOut.heightList = numpy.arange(self.processingHeaderObj.firstHeight, xf, self.processingHeaderObj.deltaHeight) + self.dataOut.channelList = range(self.systemHeaderObj.nChannels) - + self.dataOut.flagShiftFFT = True #Data is always shifted - + self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode #asumo q la data no esta decodificada - - self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip #asumo q la data esta sin flip - + + self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip #asumo q la data esta sin flip + def getData(self): """ First method to execute before "RUN" is called. - + Copia el buffer de lectura a la clase "Spectra", con todos los parametros asociados a este (metadata). cuando no hay datos en el buffer de lectura es necesario hacer una nueva lectura de los bloques de datos usando "readNextBlock" - + Return: 0 : Si no hay mas archivos disponibles 1 : Si hizo una buena copia del buffer - + Affected: self.dataOut - + self.flagDiscontinuousBlock self.flagIsNewBlock """ @@ -356,68 +367,70 @@ class SpectraReader(JRODataReader, ProcessingUnit): self.dataOut.flagNoData = True print 'Process finished' return 0 - + self.flagDiscontinuousBlock = 0 self.flagIsNewBlock = 0 - - if self.__hasNotDataInBuffer(): + + if self.__hasNotDataInBuffer(): if not( self.readNextBlock() ): self.dataOut.flagNoData = True return 0 + #data es un numpy array de 3 dmensiones (perfiles, alturas y canales) if self.data_spc is None: self.dataOut.flagNoData = True return 0 - + self.getBasicHeader() - + self.getFirstHeader() self.dataOut.data_spc = self.data_spc - + self.dataOut.data_cspc = self.data_cspc - + self.dataOut.data_dc = self.data_dc - + self.dataOut.flagNoData = False - + self.dataOut.realtime = self.online - + return self.dataOut.data_spc class SpectraWriter(JRODataWriter, Operation): - - """ + + """ Esta clase permite escribir datos de espectros a archivos procesados (.pdata). La escritura - de los datos siempre se realiza por bloques. + de los datos siempre se realiza por bloques. """ - + ext = ".pdata" - + optchar = "P" - + shape_spc_Buffer = None - + shape_cspc_Buffer = None - + shape_dc_Buffer = None - + data_spc = None - + data_cspc = None - + data_dc = None - + # dataOut = None - - def __init__(self, **kwargs): - """ + + def __init__(self): + """ Inicializador de la clase SpectraWriter para la escritura de datos de espectros. + + Affected: - Affected: self.dataOut self.basicHeaderObj self.systemHeaderObj @@ -426,49 +439,51 @@ class SpectraWriter(JRODataWriter, Operation): Return: None """ - - Operation.__init__(self, **kwargs) - + + Operation.__init__(self) + self.isConfig = False - + self.nTotalBlocks = 0 - + self.data_spc = None - + self.data_cspc = None + self.data_dc = None self.fp = None self.flagIsNewFile = 1 - - self.nTotalBlocks = 0 - + + self.nTotalBlocks = 0 + self.flagIsNewBlock = 0 self.setFile = None - + self.dtype = None - + self.path = None - + self.noMoreFiles = 0 - + self.filename = None - + self.basicHeaderObj = BasicHeader(LOCALTIME) - + self.systemHeaderObj = SystemHeader() - + self.radarControllerHeaderObj = RadarControllerHeader() - + self.processingHeaderObj = ProcessingHeader() - + def hasAllDataInBuffer(self): return 1 + def setBlockDimension(self): """ @@ -488,14 +503,15 @@ class SpectraWriter(JRODataWriter, Operation): self.shape_cspc_Buffer = (self.dataOut.nPairs, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) - + self.shape_dc_Buffer = (self.dataOut.nChannels, self.processingHeaderObj.nHeights) - + def writeBlock(self): """ Escribe el buffer en el file designado + Affected: self.data_spc @@ -504,11 +520,11 @@ class SpectraWriter(JRODataWriter, Operation): self.flagIsNewFile self.flagIsNewBlock self.nTotalBlocks - self.nWriteBlocks - + self.nWriteBlocks + Return: None """ - + spc = numpy.transpose( self.data_spc, (0,2,1) ) if not( self.processingHeaderObj.shif_fft ): spc = numpy.roll( spc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones @@ -525,6 +541,7 @@ class SpectraWriter(JRODataWriter, Operation): data['imag'] = cspc.imag data = data.reshape((-1)) data.tofile(self.fp) + if self.data_dc is not None: data = numpy.zeros( self.shape_dc_Buffer, self.dtype ) @@ -535,145 +552,147 @@ class SpectraWriter(JRODataWriter, Operation): data.tofile(self.fp) # self.data_spc.fill(0) -# +# # if self.data_dc is not None: # self.data_dc.fill(0) -# +# # if self.data_cspc is not None: # self.data_cspc.fill(0) + self.flagIsNewFile = 0 self.flagIsNewBlock = 1 self.nTotalBlocks += 1 self.nWriteBlocks += 1 self.blockIndex += 1 - + # print "[Writing] Block = %d04" %self.blockIndex - + def putData(self): """ - Setea un bloque de datos y luego los escribe en un file + Setea un bloque de datos y luego los escribe en un file + Affected: self.data_spc self.data_cspc self.data_dc - Return: - 0 : Si no hay data o no hay mas files que puedan escribirse + Return: + 0 : Si no hay data o no hay mas files que puedan escribirse 1 : Si se escribio la data de un bloque en un file """ - + if self.dataOut.flagNoData: return 0 - + self.flagIsNewBlock = 0 - + if self.dataOut.flagDiscontinuousBlock: self.data_spc.fill(0) - if self.dataOut.data_cspc is not None: - self.data_cspc.fill(0) - if self.dataOut.data_dc is not None: - self.data_dc.fill(0) + self.data_cspc.fill(0) + self.data_dc.fill(0) self.setNextFile() - + if self.flagIsNewFile == 0: self.setBasicHeader() - + self.data_spc = self.dataOut.data_spc.copy() - + if self.dataOut.data_cspc is not None: self.data_cspc = self.dataOut.data_cspc.copy() - + if self.dataOut.data_dc is not None: self.data_dc = self.dataOut.data_dc.copy() - + # #self.processingHeaderObj.dataBlocksPerFile) if self.hasAllDataInBuffer(): # self.setFirstHeader() self.writeNextBlock() - + return 1 + def __getBlockSize(self): ''' Este metodos determina el cantidad de bytes para un bloque de datos de tipo Spectra ''' - + dtype_width = self.getDtypeWidth() - + pts2write = self.dataOut.nHeights * self.dataOut.nFFTPoints - + pts2write_SelfSpectra = int(self.dataOut.nChannels * pts2write) blocksize = (pts2write_SelfSpectra*dtype_width) - + if self.dataOut.data_cspc is not None: pts2write_CrossSpectra = int(self.dataOut.nPairs * pts2write) blocksize += (pts2write_CrossSpectra*dtype_width*2) - + if self.dataOut.data_dc is not None: pts2write_DCchannels = int(self.dataOut.nChannels * self.dataOut.nHeights) blocksize += (pts2write_DCchannels*dtype_width*2) - + # blocksize = blocksize #* datatypeValue * 2 #CORREGIR ESTO return blocksize - + def setFirstHeader(self): - + """ Obtiene una copia del First Header - + Affected: self.systemHeaderObj self.radarControllerHeaderObj self.dtype - Return: + Return: None """ - + self.systemHeaderObj = self.dataOut.systemHeaderObj.copy() self.systemHeaderObj.nChannels = self.dataOut.nChannels self.radarControllerHeaderObj = self.dataOut.radarControllerHeaderObj.copy() - + self.processingHeaderObj.dtype = 1 # Spectra self.processingHeaderObj.blockSize = self.__getBlockSize() self.processingHeaderObj.profilesPerBlock = self.dataOut.nFFTPoints self.processingHeaderObj.dataBlocksPerFile = self.blocksPerFile self.processingHeaderObj.nWindows = 1 #podria ser 1 o self.dataOut.processingHeaderObj.nWindows self.processingHeaderObj.nCohInt = self.dataOut.nCohInt# Se requiere para determinar el valor de timeInterval - self.processingHeaderObj.nIncohInt = self.dataOut.nIncohInt + self.processingHeaderObj.nIncohInt = self.dataOut.nIncohInt self.processingHeaderObj.totalSpectra = self.dataOut.nPairs + self.dataOut.nChannels self.processingHeaderObj.shif_fft = self.dataOut.flagShiftFFT + if self.processingHeaderObj.totalSpectra > 0: channelList = [] for channel in range(self.dataOut.nChannels): channelList.append(channel) channelList.append(channel) - + pairsList = [] if self.dataOut.nPairs > 0: for pair in self.dataOut.pairsList: pairsList.append(pair[0]) pairsList.append(pair[1]) - + spectraComb = channelList + pairsList spectraComb = numpy.array(spectraComb, dtype="u1") self.processingHeaderObj.spectraComb = spectraComb - + if self.dataOut.code is not None: self.processingHeaderObj.code = self.dataOut.code self.processingHeaderObj.nCode = self.dataOut.nCode self.processingHeaderObj.nBaud = self.dataOut.nBaud - + if self.processingHeaderObj.nWindows != 0: self.processingHeaderObj.firstHeight = self.dataOut.heightList[0] self.processingHeaderObj.deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0] self.processingHeaderObj.nHeights = self.dataOut.nHeights self.processingHeaderObj.samplesWin = self.dataOut.nHeights - + self.processingHeaderObj.processFlags = self.getProcessFlags() - + self.setBasicHeader() diff --git a/schainpy/model/io/jroIO_voltage.py b/schainpy/model/io/jroIO_voltage.py index cc3a26f..a34c870 100644 --- a/schainpy/model/io/jroIO_voltage.py +++ b/schainpy/model/io/jroIO_voltage.py @@ -527,10 +527,16 @@ class VoltageReader(JRODataReader, ProcessingUnit): self.dataOut.flagNoData = False self.getBasicHeader() - + + #print self.basicHeaderObj.printInfo() + #print self.systemHeaderObj.printInfo() + #print self.radarControllerHeaderObj.printInfo() + #print self.processingHeaderObj.printInfo() + self.dataOut.realtime = self.online return self.dataOut.data + class VoltageWriter(JRODataWriter, Operation): """ diff --git a/schainpy/model/proc/__init__.py b/schainpy/model/proc/__init__.py index a17c2ba..af598f5 100644 --- a/schainpy/model/proc/__init__.py +++ b/schainpy/model/proc/__init__.py @@ -11,4 +11,5 @@ from jroproc_amisr import * from jroproc_correlation import * from jroproc_parameters import * from jroproc_spectra_lags import * -from jroproc_spectra_acf import * \ No newline at end of file +from jroproc_spectra_acf import * +from bltrproc_parameters import * diff --git a/schainpy/model/proc/bltrproc_parameters.py b/schainpy/model/proc/bltrproc_parameters.py new file mode 100644 index 0000000..cada61c --- /dev/null +++ b/schainpy/model/proc/bltrproc_parameters.py @@ -0,0 +1,403 @@ +''' +Created on Oct 24, 2016 + +@author: roj- LouVD +''' + +import numpy +import copy +import datetime +import time +from time import gmtime + +from numpy import transpose + +from jroproc_base import ProcessingUnit, Operation +from schainpy.model.data.jrodata import Parameters + + +class BLTRParametersProc(ProcessingUnit): + ''' + Processing unit for BLTR parameters data (winds) + + Inputs: + self.dataOut.nmodes - Number of operation modes + self.dataOut.nchannels - Number of channels + self.dataOut.nranges - Number of ranges + + self.dataOut.data_SNR - SNR array + self.dataOut.data_output - Zonal, Vertical and Meridional velocity array + self.dataOut.height - Height array (km) + self.dataOut.time - Time array (seconds) + + self.dataOut.fileIndex -Index of the file currently read + self.dataOut.lat - Latitude coordinate of BLTR location + + self.dataOut.doy - Experiment doy (number of the day in the current year) + self.dataOut.month - Experiment month + self.dataOut.day - Experiment day + self.dataOut.year - Experiment year + ''' + + def __init__(self, **kwargs): + ''' + Inputs: None + ''' + ProcessingUnit.__init__(self, **kwargs) + self.dataOut = Parameters() + self.isConfig = False + + def setup(self, mode): + ''' + ''' + self.dataOut.mode = mode + + def run(self, mode, snr_threshold=None): + ''' + Inputs: + mode = High resolution (0) or Low resolution (1) data + snr_threshold = snr filter value + ''' + + if not self.isConfig: + self.setup(mode) + self.isConfig = True + + if self.dataIn.type == 'Parameters': + self.dataOut.copy(self.dataIn) + + self.dataOut.data_output = self.dataOut.data_output[mode] + self.dataOut.heightList = self.dataOut.height[0] + self.dataOut.data_SNR = self.dataOut.data_SNR[mode] + + if snr_threshold is not None: + SNRavg = numpy.average(self.dataOut.data_SNR, axis=0) + SNRavgdB = 10*numpy.log10(SNRavg) + for i in range(3): + self.dataOut.data_output[i][SNRavgdB <= snr_threshold] = numpy.nan + +# TODO +class OutliersFilter(Operation): + + def __init__(self, **kwargs): + ''' + ''' + Operation.__init__(self, **kwargs) + + def run(self, svalue2, method, factor, filter, npoints=9): + ''' + Inputs: + svalue - string to select array velocity + svalue2 - string to choose axis filtering + method - 0 for SMOOTH or 1 for MEDIAN + factor - number used to set threshold + filter - 1 for data filtering using the standard deviation criteria else 0 + npoints - number of points for mask filter + ''' + + print ' Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor) + + + yaxis = self.dataOut.heightList + xaxis = numpy.array([[self.dataOut.utctime]]) + + # Zonal + value_temp = self.dataOut.data_output[0] + + # Zonal + value_temp = self.dataOut.data_output[1] + + # Vertical + value_temp = numpy.transpose(self.dataOut.data_output[2]) + + htemp = yaxis + std = value_temp + for h in range(len(htemp)): + nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0]) + minvalid = npoints + + #only if valid values greater than the minimum required (10%) + if nvalues_valid > minvalid: + + if method == 0: + #SMOOTH + w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1) + + + if method == 1: + #MEDIAN + w = value_temp[h] - self.Median(input=value_temp[h], width = npoints) + + dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1) + + threshold = dw*factor + value_temp[numpy.where(w > threshold),h] = numpy.nan + value_temp[numpy.where(w < -1*threshold),h] = numpy.nan + + + #At the end + if svalue2 == 'inHeight': + value_temp = numpy.transpose(value_temp) + output_array[:,m] = value_temp + + if svalue == 'zonal': + self.dataOut.data_output[0] = output_array + + elif svalue == 'meridional': + self.dataOut.data_output[1] = output_array + + elif svalue == 'vertical': + self.dataOut.data_output[2] = output_array + + return self.dataOut.data_output + + + def Median(self,input,width): + ''' + Inputs: + input - Velocity array + width - Number of points for mask filter + + ''' + + if numpy.mod(width,2) == 1: + pc = int((width - 1) / 2) + cont = 0 + output = [] + + for i in range(len(input)): + if i >= pc and i < len(input) - pc: + new2 = input[i-pc:i+pc+1] + temp = numpy.where(numpy.isfinite(new2)) + new = new2[temp] + value = numpy.median(new) + output.append(value) + + output = numpy.array(output) + output = numpy.hstack((input[0:pc],output)) + output = numpy.hstack((output,input[-pc:len(input)])) + + return output + + def Smooth(self,input,width,edge_truncate = None): + ''' + Inputs: + input - Velocity array + width - Number of points for mask filter + edge_truncate - 1 for truncate the convolution product else + + ''' + + if numpy.mod(width,2) == 0: + real_width = width + 1 + nzeros = width / 2 + else: + real_width = width + nzeros = (width - 1) / 2 + + half_width = int(real_width)/2 + length = len(input) + + gate = numpy.ones(real_width,dtype='float') + norm_of_gate = numpy.sum(gate) + + nan_process = 0 + nan_id = numpy.where(numpy.isnan(input)) + if len(nan_id[0]) > 0: + nan_process = 1 + pb = numpy.zeros(len(input)) + pb[nan_id] = 1. + input[nan_id] = 0. + + if edge_truncate == True: + output = numpy.convolve(input/norm_of_gate,gate,mode='same') + elif edge_truncate == False or edge_truncate == None: + output = numpy.convolve(input/norm_of_gate,gate,mode='valid') + output = numpy.hstack((input[0:half_width],output)) + output = numpy.hstack((output,input[len(input)-half_width:len(input)])) + + if nan_process: + pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid') + pb = numpy.hstack((numpy.zeros(half_width),pb)) + pb = numpy.hstack((pb,numpy.zeros(half_width))) + output[numpy.where(pb > 0.9999)] = numpy.nan + input[nan_id] = numpy.nan + return output + + def Average(self,aver=0,nhaver=1): + ''' + Inputs: + aver - Indicates the time period over which is averaged or consensus data + nhaver - Indicates the decimation factor in heights + + ''' + nhpoints = 48 + + lat_piura = -5.17 + lat_huancayo = -12.04 + lat_porcuya = -5.8 + + if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura: + hcm = 3. + if self.dataOut.year == 2003 : + if self.dataOut.doy >= 25 and self.dataOut.doy < 64: + nhpoints = 12 + + elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo: + hcm = 3. + if self.dataOut.year == 2003 : + if self.dataOut.doy >= 25 and self.dataOut.doy < 64: + nhpoints = 12 + + + elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya: + hcm = 5.#2 + + pdata = 0.2 + taver = [1,2,3,4,6,8,12,24] + t0 = 0 + tf = 24 + ntime =(tf-t0)/taver[aver] + ti = numpy.arange(ntime) + tf = numpy.arange(ntime) + taver[aver] + + + old_height = self.dataOut.heightList + + if nhaver > 1: + num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes + deltha = 0.05*nhaver + minhvalid = pdata*nhaver + for im in range(self.dataOut.nmodes): + new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2. + + + data_fHeigths_List = [] + data_fZonal_List = [] + data_fMeridional_List = [] + data_fVertical_List = [] + startDTList = [] + + + for i in range(ntime): + height = old_height + + start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5) + stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5) + + + limit_sec1 = time.mktime(start.timetuple()) + limit_sec2 = time.mktime(stop.timetuple()) + + t1 = numpy.where(self.f_timesec >= limit_sec1) + t2 = numpy.where(self.f_timesec < limit_sec2) + time_select = [] + for val_sec in t1[0]: + if val_sec in t2[0]: + time_select.append(val_sec) + + + time_select = numpy.array(time_select,dtype = 'int') + minvalid = numpy.ceil(pdata*nhpoints) + + zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan + mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan + ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan + + if nhaver > 1: + new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan + new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan + new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan + + if len(time_select) > minvalid: + time_average = self.f_timesec[time_select] + + for im in range(self.dataOut.nmodes): + + for ih in range(self.dataOut.nranges): + if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid: + zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) + + if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid: + mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) + + if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid: + ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) + + if nhaver > 1: + for ih in range(num_hei): + hvalid = numpy.arange(nhaver) + nhaver*ih + + if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid: + new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) / numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) + + if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid: + new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) / numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) + + if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid: + new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) / numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) + if nhaver > 1: + zon_aver = new_zon_aver + mer_aver = new_mer_aver + ver_aver = new_ver_aver + height = new_height + + + tstart = time_average[0] + tend = time_average[-1] + startTime = time.gmtime(tstart) + + year = startTime.tm_year + month = startTime.tm_mon + day = startTime.tm_mday + hour = startTime.tm_hour + minute = startTime.tm_min + second = startTime.tm_sec + + startDTList.append(datetime.datetime(year,month,day,hour,minute,second)) + + + o_height = numpy.array([]) + o_zon_aver = numpy.array([]) + o_mer_aver = numpy.array([]) + o_ver_aver = numpy.array([]) + if self.dataOut.nmodes > 1: + for im in range(self.dataOut.nmodes): + + if im == 0: + h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:]))) + else: + h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:]))) + + + ht = h_select[0] + + o_height = numpy.hstack((o_height,height[im,ht])) + o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im])) + o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im])) + o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im])) + + data_fHeigths_List.append(o_height) + data_fZonal_List.append(o_zon_aver) + data_fMeridional_List.append(o_mer_aver) + data_fVertical_List.append(o_ver_aver) + + + else: + h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:]))) + ht = h_select[0] + o_height = numpy.hstack((o_height,height[im,ht])) + o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im])) + o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im])) + o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im])) + + data_fHeigths_List.append(o_height) + data_fZonal_List.append(o_zon_aver) + data_fMeridional_List.append(o_mer_aver) + data_fVertical_List.append(o_ver_aver) + + + return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List + + + \ No newline at end of file diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 299b3d4..5df7042 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -901,3 +901,4 @@ class IncohInt(Operation): dataOut.nIncohInt *= self.n dataOut.utctime = avgdatatime dataOut.flagNoData = False + diff --git a/schainpy/model/utils/jroutils_publish.py b/schainpy/model/utils/jroutils_publish.py index 726065a..cea868b 100644 --- a/schainpy/model/utils/jroutils_publish.py +++ b/schainpy/model/utils/jroutils_publish.py @@ -113,6 +113,8 @@ class Data(object): self.__heights = [] self.__all_heights = set() for plot in self.plottypes: + if 'snr' in plot: + plot = 'snr' self.data[plot] = {} def shape(self, key): @@ -138,8 +140,9 @@ class Data(object): self.parameters = getattr(dataOut, 'parameters', []) self.pairs = dataOut.pairsList self.channels = dataOut.channelList - self.xrange = (dataOut.getFreqRange(1)/1000. , dataOut.getAcfRange(1) , dataOut.getVelRange(1)) self.interval = dataOut.getTimeInterval() + if 'spc' in self.plottypes or 'cspc' in self.plottypes: + self.xrange = (dataOut.getFreqRange(1)/1000. , dataOut.getAcfRange(1) , dataOut.getVelRange(1)) self.__heights.append(dataOut.heightList) self.__all_heights.update(dataOut.heightList) self.__times.append(tm) diff --git a/schainpy/project.py b/schainpy/project.py deleted file mode 100644 index bdea4c3..0000000 --- a/schainpy/project.py +++ /dev/null @@ -1,34 +0,0 @@ -from schainpy.controller import Project - -desc = "A schain project" - -controller = Project() -controller.setup(id='191', name="project", description=desc) - -readUnitConf = controller.addReadUnit(datatype='VoltageReader', - path="/home/nanosat/schain/schainpy", - startDate="1970/01/01", - endDate="2017/12/31", - startTime="00:00:00", - endTime="23:59:59", - online=0, - verbose=1, - walk=1, - ) - -procUnitConf1 = controller.addProcUnit(datatype='VoltageProc', inputId=readUnitConf.getId()) - -opObj11 = procUnitConf1.addOperation(name='ProfileSelector', optype='other') -opObj11.addParameter(name='profileRangeList', value='120,183', format='intlist') - -opObj11 = procUnitConf1.addOperation(name='RTIPlot', optype='other') -opObj11.addParameter(name='wintitle', value='Jicamarca Radio Observatory', format='str') -opObj11.addParameter(name='showprofile', value='0', format='int') -opObj11.addParameter(name='xmin', value='0', format='int') -opObj11.addParameter(name='xmax', value='24', format='int') -opObj11.addParameter(name='figpath', value="/home/nanosat/schain/schainpy/figs", format='str') -opObj11.addParameter(name='wr_period', value='5', format='int') -opObj11.addParameter(name='exp_code', value='22', format='int') - - -controller.start() diff --git a/schainpy/utils/trash b/schainpy/utils/trash deleted file mode 100644 index 384299d..0000000 --- a/schainpy/utils/trash +++ /dev/null @@ -1 +0,0 @@ -You should install "digital_rf_hdf5" module if you want to read USRP data