From 3cc55c17931803216e5276b27591bd9753f9f6b1 2015-08-06 17:47:33 From: Julio Valdez Date: 2015-08-06 17:47:33 Subject: [PATCH] data -Corrections in Parameters class graphics -Bug fix in RTI Plot -Bug fix in Sky Map Plot -New plot added: Phase Plot -New plotting codes added for new graphics model -Bug fix in HDF5 Writing module -Added new operation and plot to calculate Channel Phase Offset in JASMET experiments scripts -New scripts for JASMET campaign --- diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 5c44d60..0764574 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -1069,13 +1069,13 @@ class Parameters(JROData): data_SNR = None #Signal to Noise Ratio - heightRange = None #Heights +# heightRange = None #Heights - abscissaRange = None #Abscissa, can be velocities, lags or time + abscissaList = None #Abscissa, can be velocities, lags or time noise = None #Noise Potency -# initUtcTime = None #Initial UTC time + utctimeInit = None #Initial UTC time paramInterval = None #Time interval to calculate Parameters in seconds @@ -1093,6 +1093,8 @@ class Parameters(JROData): data_output = None #Out signal + + def __init__(self): ''' Constructor @@ -1102,14 +1104,21 @@ class Parameters(JROData): self.systemHeaderObj = SystemHeader() self.type = "Parameters" - + def getTimeRange1(self): datatime = [] - datatime.append(self.ltctime) - datatime.append(self.ltctime + self.outputInterval - 1) + if self.useLocalTime: + time1 = self.utctimeInit - self.timeZone*60 + else: + time1 = utctimeInit + +# datatime.append(self.utctimeInit) +# datatime.append(self.utctimeInit + self.outputInterval) + datatime.append(time1) + datatime.append(time1 + self.outputInterval) datatime = numpy.array(datatime) - return datatime + return datatime diff --git a/schainpy/model/graphics/figure.py b/schainpy/model/graphics/figure.py index 2b3f63b..7369adb 100644 --- a/schainpy/model/graphics/figure.py +++ b/schainpy/model/graphics/figure.py @@ -84,9 +84,9 @@ class Figure(Operation): #raise ValueError, "(timerange) or (xmin & xmax) should be defined" if timerange != None: - txmin = x[0] - x[0] % min(timerange/10, 10*60) + txmin = x[0] #- x[0] % min(timerange/10, 10*60) else: - txmin = x[0] - x[0] % 10*60 + txmin = x[0] #- x[0] % 10*60 thisdatetime = datetime.datetime.utcfromtimestamp(txmin) thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0)) diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index 260e844..de8969c 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -198,7 +198,7 @@ class SkyMapPlot(Figure): WIDTHPROF = None HEIGHTPROF = None - PREFIX = 'prm' + PREFIX = 'mmap' def __init__(self): @@ -213,7 +213,7 @@ class SkyMapPlot(Figure): self.HEIGHTPROF = 0 self.counter_imagwr = 0 - self.PLOT_CODE = SKYMAP_CODE + self.PLOT_CODE = MSKYMAP_CODE self.FTP_WEI = None self.EXP_CODE = None @@ -269,7 +269,7 @@ class SkyMapPlot(Figure): zmax : None """ - arrayParameters = dataOut.data_param + arrayParameters = dataOut.data_param[0,:] error = arrayParameters[:,-1] indValid = numpy.where(error == 0)[0] finalMeteor = arrayParameters[indValid,:] @@ -324,6 +324,7 @@ class SkyMapPlot(Figure): ftp=ftp, wr_period=wr_period, thisDatetime=thisDatetime) + class WindProfilerPlot(Figure): @@ -336,7 +337,7 @@ class WindProfilerPlot(Figure): def __init__(self): - self.timerange = 2*60*60 + self.timerange = None self.__isConfig = False self.__nsubplots = 1 @@ -391,7 +392,7 @@ class WindProfilerPlot(Figure): self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1) counter += 1 - def run(self, dataOut, id, wintitle="", channelList=None, + def run(self, dataOut, id, wintitle="", channelList=None, showprofile='False', xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None, timerange=None, SNRthresh = None, @@ -423,11 +424,11 @@ class WindProfilerPlot(Figure): raise ValueError, "Channel %d is not in dataOut.channelList" channelIndexList.append(dataOut.channelList.index(channel)) - if timerange != None: - self.timerange = timerange - - tmin = None - tmax = None +# if timerange != None: +# self.timerange = timerange +# +# tmin = None +# tmax = None x = dataOut.getTimeRange1() # y = dataOut.heightList @@ -453,7 +454,7 @@ class WindProfilerPlot(Figure): z[i,ind] = numpy.nan - showprofile = False +# showprofile = False # thisDatetime = dataOut.datatime thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + "Wind" @@ -468,6 +469,9 @@ class WindProfilerPlot(Figure): showprofile=showprofile, show=show) + if timerange != None: + self.timerange = timerange + self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange) if ymin == None: ymin = numpy.nanmin(y) @@ -485,17 +489,18 @@ class WindProfilerPlot(Figure): if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB) if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB) + self.FTP_WEI = ftp_wei self.EXP_CODE = exp_code self.SUB_EXP_CODE = sub_exp_code self.PLOT_POS = plot_pos - + self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") self.__isConfig = True - - + self.figfile = figfile + self.setWinTitle(title) - + if ((self.xmax - x[1]) < (x[1]-x[0])): x[1] = self.xmax @@ -620,7 +625,7 @@ class ParametersPlot(Figure): def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False, xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None, parameterIndex = None, onlyPositive = False, - SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, + SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, onlySNR = False, DOP = True, zlabel = "", parameterName = "", parameterObject = "data_param", save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True, @@ -689,12 +694,14 @@ class ParametersPlot(Figure): # SNRavgdB = 10*numpy.log10(SNRavg) ind = numpy.where(SNRdB < 10**(SNRthresh/10)) z[ind] = numpy.nan - + 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 + if onlyPositive: colormap = "jet" zmin = 0 @@ -734,9 +741,23 @@ class ParametersPlot(Figure): x[1] = self.xmax for i in range(nchan): + + if (SNR and not onlySNR): j = 2*i + else: j = i j = nGraphsByChannel*i + + if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)): + title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith) + if not onlySNR: + axes = self.axesList[j*self.__nsubplots] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap, + ticksize=9, cblabel=zlabel, cbsize="1%") + if DOP: title = "%s Channel %d: %s" %(parameterName, channelIndexList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) @@ -751,7 +772,12 @@ class ParametersPlot(Figure): if SNR: title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i]+1, 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)] + 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, @@ -1155,10 +1181,182 @@ class EWDriftsPlot(Figure): self.__isConfig = False self.figfile = None + + + +class PhasePlot(Figure): + + __isConfig = None + __nsubplots = None + + PREFIX = 'mphase' + + def __init__(self): + + self.timerange = 24*60*60 + self.__isConfig = False + self.__nsubplots = 1 + self.counter_imagwr = 0 + self.WIDTH = 600 + self.HEIGHT = 300 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.xdata = None + self.ydata = None + + self.PLOT_CODE = MPHASE_CODE + + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + + + self.filename_phase = None + + self.figfile = None + + 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="", pairsList=None, showprofile='True', + xmin=None, xmax=None, ymin=None, ymax=None, + timerange=None, + save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): + + + if timerange != None: + self.timerange = timerange + + tmin = None + tmax = None + x = dataOut.getTimeRange1() + y = dataOut.getHeiRange() + + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y")) + xlabel = "Local Time" + ylabel = "Phase" + + + #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList))) + phase_beacon = dataOut.data_output + + + if not self.__isConfig: + + self.nplots = phase_beacon.size + + self.setup(id=id, + nplots=self.nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + tmin, tmax = self.getTimeLim(x, xmin, xmax) + if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0 + if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0 + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") + self.__isConfig = True + self.figfile = figfile + self.xdata = numpy.array([]) + self.ydata = numpy.array([]) + + #open file beacon phase + path = '%s%03d' %(self.PREFIX, self.id) + beacon_file = os.path.join(path,'%s.txt'%self.name) + self.filename_phase = os.path.join(figpath,beacon_file) + #self.save_phase(self.filename_phase) + + + #store data beacon phase + #self.save_data(self.filename_phase, phase_beacon, thisDatetime) + + self.setWinTitle(title) + + + title = "Phase Offset %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + + legendlabels = ["phase %d"%(chan) for chan in numpy.arange(self.nplots)] + + axes = self.axesList[0] + + self.xdata = numpy.hstack((self.xdata, x[0:1])) + + if len(self.ydata)==0: + self.ydata = phase_beacon.reshape(-1,1) + else: + self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1))) + + + axes.pmultilineyaxis(x=self.xdata, y=self.ydata, + xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid", + XAxisAsTime=True, grid='both' + ) + + self.draw() + + if x[1] >= self.axesList[0].xmax: + self.counter_imagwr = wr_period + del self.xdata + del self.ydata + self.__isConfig = False + + if self.figfile == None: + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + self.figfile = self.getFilename(name = str_datetime) + + if figpath != '': + self.counter_imagwr += 1 + if (self.counter_imagwr>=wr_period): + # store png plot to local folder + self.saveFigure(figpath, self.figfile) + # store png plot to FTP server according to RT-Web format + name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) + ftp_filename = os.path.join(figpath, name) + self.saveFigure(figpath, ftp_filename) + self.counter_imagwr = 0 + self.figfile = None + self.save(figpath=figpath, figfile=figfile, save=save, ftp=ftp, wr_period=wr_period, thisDatetime=thisDatetime, - update_figfile=False) \ No newline at end of file + update_figfile=False) diff --git a/schainpy/model/graphics/mpldriver.py b/schainpy/model/graphics/mpldriver.py index 3091b6a..676791d 100644 --- a/schainpy/model/graphics/mpldriver.py +++ b/schainpy/model/graphics/mpldriver.py @@ -395,9 +395,11 @@ def createPolar(ax, x, y, # 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, -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') - printLabels(ax, xlabel, '', title) + ax.yaxis.labelpad = 230 + printLabels(ax, xlabel, ylabel, title) iplot = ax.lines[-1] if '0.' in matplotlib.__version__[0:2]: @@ -423,7 +425,7 @@ 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, '', title) + printLabels(ax, xlabel, ylabel, title) set_linedata(ax, x, y, idline=0) diff --git a/schainpy/model/graphics/plotting_codes.py b/schainpy/model/graphics/plotting_codes.py index 5fc8fa8..3432f0c 100644 --- a/schainpy/model/graphics/plotting_codes.py +++ b/schainpy/model/graphics/plotting_codes.py @@ -18,10 +18,11 @@ NOISE_CODE = 17 BEACON_CODE = 18 #USED IN jroplot_parameters.py - -MOMENTS_CODE = 20 -SKYMAP_CODE = 21 WIND_CODE = 22 -PARMS_CODE = 23 -SPECFIT_CODE = 24 -EWDRIFT_CODE = 25 +MSKYMAP_CODE = 23 +MPHASE_CODE = 24 + +MOMENTS_CODE = 25 +PARMS_CODE = 26 +SPECFIT_CODE = 27 +EWDRIFT_CODE = 28 diff --git a/schainpy/model/io/jroIO_HDF5.py b/schainpy/model/io/jroIO_HDF5.py index c3c4fa2..50ae899 100644 --- a/schainpy/model/io/jroIO_HDF5.py +++ b/schainpy/model/io/jroIO_HDF5.py @@ -3,6 +3,7 @@ import time import os import h5py import re +import tables from schainpy.model.data.jrodata import * from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation @@ -156,7 +157,7 @@ class HDF5Reader(ProcessingUnit): for thisPath in os.listdir(single_path): if not os.path.isdir(os.path.join(single_path,thisPath)): continue - if not isRadarFolder(thisPath): + if not isDoyFolder(thisPath): continue dirList.append(thisPath) @@ -549,6 +550,8 @@ class HDF5Writer(Operation): metaFile = None + filename = None + path = None setFile = None @@ -586,9 +589,11 @@ class HDF5Writer(Operation): mode = None nDatas = None #Number of datasets to be stored per array - + nDims = None #Number Dimensions in each dataset + nDimsForDs = None + def __init__(self): Operation.__init__(self) @@ -675,16 +680,26 @@ class HDF5Writer(Operation): setFile = self.setFile timeTuple = time.localtime(self.dataOut.utctime) - subfolder = '' - + + subfolder = '' fullpath = os.path.join( path, subfolder ) + if not( os.path.exists(fullpath) ): os.mkdir(fullpath) setFile = -1 #inicializo mi contador de seteo + + subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday) + fullpath = os.path.join( path, subfolder ) + + if not( os.path.exists(fullpath) ): + os.mkdir(fullpath) + setFile = -1 #inicializo mi contador de seteo + else: filesList = os.listdir( fullpath ) + filesList = sorted( filesList, key=str.lower ) if len( filesList ) > 0: - filesList = sorted( filesList, key=str.lower ) + filesList = [k for k in filesList if 'M' in k] filen = filesList[-1] # el filename debera tener el siguiente formato # 0 1234 567 89A BCDE (hex) @@ -697,7 +712,7 @@ class HDF5Writer(Operation): setFile = -1 #inicializo mi contador de seteo setFile += 1 - + file = '%s%4.4d%3.3d%3.3d%s' % (self.metaoptchar, timeTuple.tm_year, timeTuple.tm_yday, @@ -726,19 +741,15 @@ class HDF5Writer(Operation): path = self.path setFile = self.setFile mode = self.mode - - if self.fp != None: - self.fp.close() - + timeTuple = time.localtime(self.dataOut.utctime) subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday) fullpath = os.path.join( path, subfolder ) - if not( os.path.exists(fullpath) ): - os.mkdir(fullpath) - setFile = -1 #inicializo mi contador de seteo - else: + + if os.path.exists(fullpath): filesList = os.listdir( fullpath ) + filesList = [k for k in filesList if 'D' in k] if len( filesList ) > 0: filesList = sorted( filesList, key=str.lower ) filen = filesList[-1] @@ -771,17 +782,22 @@ class HDF5Writer(Operation): ds = [] data = [] + nDimsForDs = [] nDatas = numpy.zeros(len(self.dataList)) nDims = self.arrayDim[:,0] + nDim1 = self.arrayDim[:,2] + nDim0 = self.arrayDim[:,3] + for i in range(len(self.dataList)): if nDims[i]==1: - ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,None) , chunks = True, dtype='S20') +# ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype='S20') + ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype=numpy.float64) ds.append(ds0) data.append([]) - + nDimsForDs.append(nDims[i]) else: if mode[i]==0: @@ -800,20 +816,22 @@ class HDF5Writer(Operation): tableName = strMode + str(j) if nDims[i] == 3: - ds0 = grp0.create_dataset(tableName, (1,1,1) , maxshape=(None,None,None), chunks=True) + ds0 = grp0.create_dataset(tableName, (nDim1[i],nDim0[i],1) , data = numpy.zeros((nDim1[i],nDim0[i],1)) ,maxshape=(None,nDim0[i],None), chunks=True) else: - ds0 = grp0.create_dataset(tableName, (1,1) , maxshape=(None,None), chunks=True) + ds0 = grp0.create_dataset(tableName, (1,nDim0[i]), data = numpy.zeros((1,nDim0[i])) , maxshape=(None,nDim0[i]), chunks=True) ds.append(ds0) data.append([]) - + nDimsForDs.append(nDims[i]) self.nDatas = nDatas self.nDims = nDims - + self.nDimsForDs = nDimsForDs #Saving variables print 'Writing the file: %s'%filename + self.filename = filename self.fp = fp self.grp = grp + self.grp.attrs.modify('nRecords', 1) self.ds = ds self.data = data @@ -823,13 +841,65 @@ class HDF5Writer(Operation): return def putData(self): + + if not self.firsttime: + self.fp.flush() + self.fp.close() + self.readBlock() + + if self.blockIndex == self.blocksPerFile: + + self.setNextFile() + self.setBlock() self.writeBlock() - if self.blockIndex == self.blocksPerFile: - self.setNextFile() + return + def readBlock(self): + + ''' + data Array configured + + + self.data + ''' + ds = self.ds + #Setting HDF5 File + fp = h5py.File(self.filename,'r+') + grp = fp["Data"] + ind = 0 + +# grp.attrs['blocksPerFile'] = 0 + for i in range(len(self.dataList)): + + if self.nDims[i]==1: + ds0 = grp[self.dataList[i]] + ds[ind] = ds0 + ind += 1 + else: + if self.mode[i]==0: + strMode = "channel" + else: + strMode = "param" + + grp0 = grp[self.dataList[i]] + + for j in range(int(self.nDatas[i])): + tableName = strMode + str(j) + ds0 = grp0[tableName] + ds[ind] = ds0 + ind += 1 + + + self.fp = fp + self.grp = grp + self.ds = ds + + return + + def setBlock(self): ''' data Array configured @@ -848,11 +918,11 @@ class HDF5Writer(Operation): dataAux = getattr(self.dataOut,self.dataList[i]) if nDims[i] == 1: - data[ind] = numpy.array([str(dataAux)]).reshape((1,1)) - if not self.firsttime: - data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind])) +# data[ind] = numpy.array([str(dataAux)]).reshape((1,1)) + data[ind] = dataAux +# if not self.firsttime: +# data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind])) ind += 1 - else: for j in range(int(nDatas[i])): if (mode[i] == 0) or (nDims[i] == 2): #In case division per channel or Dimensions is only 1 @@ -860,19 +930,19 @@ class HDF5Writer(Operation): else: data[ind] = dataAux[:,j,:] - if nDims[i] == 3: - data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1)) +# if nDims[i] == 3: +# data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1)) - if not self.firsttime: - data[ind] = numpy.dstack((self.ds[ind][:], data[ind])) +# if not self.firsttime: +# data[ind] = numpy.dstack((self.ds[ind][:], data[ind])) - else: - data[ind] = data[ind].reshape((1,data[ind].shape[0])) +# else: +# data[ind] = data[ind].reshape((1,data[ind].shape[0])) - if not self.firsttime: - data[ind] = numpy.vstack((self.ds[ind][:], data[ind])) +# if not self.firsttime: +# data[ind] = numpy.vstack((self.ds[ind][:], data[ind])) ind += 1 - + self.data = data return @@ -881,13 +951,45 @@ class HDF5Writer(Operation): Saves the block in the HDF5 file ''' for i in range(len(self.ds)): - self.ds[i].resize(self.data[i].shape) - self.ds[i][:] = self.data[i] + if self.firsttime: +# self.ds[i].resize(self.data[i].shape) +# self.ds[i][self.blockIndex,:] = self.data[i] + if type(self.data[i]) == numpy.ndarray: + nDims1 = len(self.ds[i].shape) + + if nDims1 == 3: + self.data[i] = self.data[i].reshape((self.data[i].shape[0],self.data[i].shape[1],1)) + + self.ds[i].resize(self.data[i].shape) + self.ds[i][:] = self.data[i] + else: + if self.nDimsForDs[i] == 1: + self.ds[i].resize((self.ds[i].shape[0], self.ds[i].shape[1] + 1)) + self.ds[i][0,-1] = self.data[i] + elif self.nDimsForDs[i] == 2: + self.ds[i].resize((self.ds[i].shape[0] + 1,self.ds[i].shape[1])) + self.ds[i][self.blockIndex,:] = self.data[i] + elif self.nDimsForDs[i] == 3: + + dataShape = self.data[i].shape + dsShape = self.ds[i].shape + + if dataShape[0]==dsShape[0]: + self.ds[i].resize((self.ds[i].shape[0],self.ds[i].shape[1],self.ds[i].shape[2]+1)) + self.ds[i][:,:,-1] = self.data[i] + else: + self.ds[i].resize((self.ds[i].shape[0] + dataShape[0],self.ds[i].shape[1],self.ds[i].shape[2])) + self.ds[i][dsShape[0]:,:,0] = self.data[i] +# self.ds[i].append(self.data[i]) +# self.fp.flush() +# if not self.firsttime: +# self.fp.root.Data._v_attrs.nRecords = self.blockIndex + +# if self.firsttime: +# self.fp.close() +# self.readBlock2() self.blockIndex += 1 - - self.grp.attrs.modify('nRecords', self.blockIndex) - self.firsttime = False return diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index b64eb73..1c6c1fa 100644 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -239,13 +239,34 @@ class ParametersProc(ProcessingUnit): return moments #------------------ Get SA Parameters -------------------------- - def GetSAParameters(self): - data = self.dataOut.data_pre - crossdata = self.dataIn.data_cspc - a = 1 - - + def GetSAParameters(self): + pairslist = self.dataOut.groupList + num_pairs = len(pairslist) + + vel = self.dataOut.abscissaList + spectra = self.dataOut.data_pre + cspectra = self.dataIn.data_cspc + delta_v = vel[1] - vel[0] + + #Calculating the power spectrum + spc_pow = numpy.sum(spectra, 3)*delta_v + #Normalizing Spectra + norm_spectra = spectra/spc_pow + #Calculating the norm_spectra at peak + max_spectra = numpy.max(norm_spectra, 3) + + #Normalizing Cross Spectra + norm_cspectra = numpy.zeros(cspectra.shape) + + for i in range(num_chan): + norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:]) + + max_cspectra = numpy.max(norm_cspectra,2) + max_cspectra_index = numpy.argmax(norm_cspectra, 2) + + for i in range(num_pairs): + cspc_par[i,:,:] = __calculateMoments(norm_cspectra) #------------------- Get Lags ---------------------------------- def GetLags(self): @@ -339,8 +360,8 @@ class ParametersProc(ProcessingUnit): return phase #------------------- Detect Meteors ------------------------------ - def DetectMeteors(self, hei_ref = None, tauindex = 0, - predefinedPhaseShifts = None, centerReceiverIndex = 2, + def MeteorDetection(self, hei_ref = None, tauindex = 0, + predefinedPhaseShifts = None, centerReceiverIndex = 2, saveAll = False, cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, noise_timeStep = 4, noise_multiple = 4, multDet_timeLimit = 1, multDet_rangeLimit = 3, @@ -430,14 +451,21 @@ class ParametersProc(ProcessingUnit): if predefinedPhaseShifts != None: hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180 - else: + + elif beaconPhaseShifts: #get hardware phase shifts using beacon signal hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10) hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0) - + + else: + hardwarePhaseShifts = numpy.zeros(5) + + voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex') for i in range(self.dataOut.data_pre.shape[0]): voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i]) + + #******************END OF REMOVING HARDWARE PHASE DIFFERENCES********* #Remove DC @@ -497,29 +525,40 @@ class ParametersProc(ProcessingUnit): listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval) if len(listMeteors4) > 0: - #Setting New Array - date = repr(self.dataOut.datatime) - arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) - #Calculate AOA (Error N 3, 4) - #JONES ET AL. 1998 - AOAthresh = numpy.pi/8 - error = arrayParameters[:,-1] - phases = -arrayMeteors4[:,9:13] pairsList = [] - pairsList.append((0,3)) - pairsList.append((1,2)) - arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth) - - #Calculate Heights (Error N 13 and 14) - error = arrayParameters[:,-1] - Ranges = arrayParameters[:,2] - zenith = arrayParameters[:,5] - arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax) + pairx = (0,3) + pairy = (1,2) + pairsList.append(pairx) + pairsList.append(pairy) + + #Setting New Array + date = repr(self.dataOut.datatime) + arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) + + meteorOps = MeteorOperations() + jph = numpy.array([0,0,0,0]) + h = (hmin,hmax) + arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph) + +# #Calculate AOA (Error N 3, 4) +# #JONES ET AL. 1998 +# error = arrayParameters[:,-1] +# AOAthresh = numpy.pi/8 +# phases = -arrayParameters[:,9:13] +# arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth) +# +# #Calculate Heights (Error N 13 and 14) +# error = arrayParameters[:,-1] +# Ranges = arrayParameters[:,2] +# zenith = arrayParameters[:,5] +# arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax) +# error = arrayParameters[:,-1] #********************* END OF PARAMETERS CALCULATION ************************** - #***************************+ SAVE DATA IN HDF5 FORMAT ********************** - self.dataOut.data_param = arrayParameters + #***************************+ PASS DATA TO NEXT STEP ********************** + arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1])) + self.dataOut.data_param = arrayFinal return @@ -948,9 +987,8 @@ class ParametersProc(ProcessingUnit): listMeteors1 = [] for i in range(len(listMeteors)): - meteor = listMeteors[i] - meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1])) - if meteor[-1] == 0: + meteorAux = listMeteors[i] + if meteorAux[-1] == 0: mStart = listMeteors[i][1] mPeak = listMeteors[i][2] mLag = mPeak - mStart + lag @@ -1007,7 +1045,7 @@ class ParametersProc(ProcessingUnit): #New arrays arrayMeteors = numpy.array(listMeteors) - arrayParameters = numpy.zeros((len(listMeteors),10)) + arrayParameters = numpy.zeros((len(listMeteors), 14)) #Date inclusion date = re.findall(r'\((.*?)\)', date) @@ -1017,14 +1055,18 @@ class ParametersProc(ProcessingUnit): arrayDate = numpy.tile(date, (len(listMeteors), 1)) #Meteor array - arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] - arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) +# arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] +# arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) #Parameters Array - arrayParameters[:,0:3] = arrayMeteors[:,0:3] - arrayParameters[:,-3:] = arrayMeteors[:,-3:] - - return arrayMeteors, arrayParameters + arrayParameters[:,:2] = arrayDate #Date + arrayParameters[:,2] = heiRang[arrayMeteors[:,0].astype(int)] #Range + arrayParameters[:,7:9] = arrayMeteors[:,-3:-1] #Radial velocity and its error + arrayParameters[:,9:13] = arrayMeteors[:,7:11] #Phases + arrayParameters[:,-1] = arrayMeteors[:,-1] #Error + + + return arrayParameters def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth): @@ -1246,7 +1288,6 @@ class ParametersProc(ProcessingUnit): self.dataOut.data_param[i,:,h] = minp return - def __residFunction(self, p, dp, LT, constants): fm = self.dataOut.library.modelFunction(p, constants) @@ -1327,7 +1368,6 @@ class WindProfiler(Operation): return A1 def __correctValues(self, heiRang, phi, velRadial, SNR): - listPhi = phi.tolist() maxid = listPhi.index(max(listPhi)) minid = listPhi.index(min(listPhi)) @@ -1539,14 +1579,15 @@ class WindProfiler(Operation): Output: Winds estimation (Zonal and Meridional) Parameters affected: Winds - ''' + ''' + print arrayMeteor.shape #Settings nInt = (heightMax - heightMin)/2 winds = numpy.zeros((2,nInt))*numpy.nan #Filter errors - error = numpy.where(arrayMeteor[:,-1] == 0)[0] - finalMeteor = arrayMeteor[error,:] + error = numpy.where(arrayMeteor[0,:,-1] == 0)[0] + finalMeteor = arrayMeteor[0,error,:] #Meteor Histogram finalHeights = finalMeteor[:,3] @@ -1592,10 +1633,10 @@ class WindProfiler(Operation): def run(self, dataOut, technique, **kwargs): param = dataOut.data_param - if dataOut.abscissaList != None: - absc = dataOut.abscissaList[:-1] +# if dataOut.abscissaList != None: +# absc = dataOut.abscissaList[:-1] noise = dataOut.noise - heightList = dataOut.getHeiRange() + heightList = dataOut.heightList SNR = dataOut.data_SNR if technique == 'DBS': @@ -1616,19 +1657,14 @@ class WindProfiler(Operation): else: correctFactor = 1 if kwargs.has_key('channelList'): channelList = kwargs['channelList'] + if len(channelList) == 2: + horizontalOnly = True arrayChannel = numpy.array(channelList) param = param[arrayChannel,:,:] theta_x = theta_x[arrayChannel] theta_y = theta_y[arrayChannel] velRadial0 = param[:,1,:] #Radial velocity - - if velRadial0.shape[0] != theta_x.shape[0] or velRadial0.shape[0] != theta_y.shape[0]: - raise ValueError, "The max number of channels is %d, and the length of cosine director is %d. Please check: dirCosX, dirCosY, elevation or azimuth arguments" %(velRadial0.shape[0], theta_x.shape[0]) - - if theta_x.shape[0] == 2: - horizontalOnly = True - dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.timeInterval @@ -1685,7 +1721,7 @@ class WindProfiler(Operation): if self.__isConfig == False: # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) #Get Initial LTC time - self.__initime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime) + self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() self.__isConfig = True @@ -1695,7 +1731,7 @@ class WindProfiler(Operation): self.__firstdata = copy.copy(dataOut) else: - self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) + self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param)) self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready @@ -1779,8 +1815,330 @@ class EWDriftsEstimation(Operation): dataOut.outputInterval = dataOut.timeInterval return +class PhaseCalibration(Operation): + + __buffer = None + + __initime = None + __dataReady = False + + __isConfig = False + + def __checkTime(self, currentTime, initTime, paramInterval, outputInterval): + + dataTime = currentTime + paramInterval + deltaTime = dataTime - initTime + + if deltaTime >= outputInterval or deltaTime < 0: + return True + + return False + + def __getGammas(self, pairs, k, d, phases): + gammas = numpy.zeros(2) + + for i in range(len(pairs)): + + pairi = pairs[i] + + #Calculating gamma + jdcos = phases[:,pairi[1]]/(k*d[pairi[1]]) + jgamma = numpy.angle(numpy.exp(1j*(k*d[pairi[0]]*jdcos - phases[:,pairi[0]]))) + + #Revised distribution + jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi)) + + #Histogram + nBins = 64.0 + rmin = -0.5*numpy.pi + rmax = 0.5*numpy.pi + phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax)) + + meteorsY = phaseHisto[0] + phasesX = phaseHisto[1][:-1] + width = phasesX[1] - phasesX[0] + phasesX += width/2 + + #Gaussian aproximation + bpeak = meteorsY.argmax() + peak = meteorsY.max() + jmin = bpeak - 5 + jmax = bpeak + 5 + 1 + + if jmin<0: + jmin = 0 + jmax = 6 + elif jmax > meteorsY.size: + jmin = meteorsY.size - 6 + jmax = meteorsY.size + + x0 = numpy.array([peak,bpeak,50]) + coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax])) + + #Gammas + gammas[i] = coeff[0][1] +# gammas[i] = bpeak + + return gammas + + def __residualFunction(self, coeffs, y, t): + + return y - self.__gauss_function(t, coeffs) + + def __gauss_function(self, t, coeffs): + + return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2) + + def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray): + meteorOps = MeteorOperations() + nchan = 4 + pairx = pairsList[0] + pairy = pairsList[1] + center_xangle = 0 + center_yangle = 0 + range_angle = numpy.array([8*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4]) + ntimes = len(range_angle) + + nstepsx = 20.0 + nstepsy = 20.0 + + for iz in range(ntimes): + min_xangle = -range_angle[iz]/2 + center_xangle + max_xangle = range_angle[iz]/2 + center_xangle + min_yangle = -range_angle[iz]/2 + center_yangle + max_yangle = range_angle[iz]/2 + center_yangle + + inc_x = (max_xangle-min_xangle)/nstepsx + inc_y = (max_yangle-min_yangle)/nstepsy + + alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle + alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle + penalty = numpy.zeros((nstepsx,nstepsy)) + jph_array = numpy.zeros((nchan,nstepsx,nstepsy)) + jph = numpy.zeros(nchan) + + # Iterations looking for the offset + for iy in range(int(nstepsy)): + for ix in range(int(nstepsx)): + jph[pairy[1]] = alpha_y[iy] + jph[pairy[0]] = -gammas[1] + alpha_y[iy]*d[pairy[0]]/d[pairy[1]] + + jph[pairx[1]] = alpha_x[ix] + jph[pairx[0]] = -gammas[0] + alpha_x[ix]*d[pairx[0]]/d[pairx[1]] + + jph_array[:,ix,iy] = jph + + meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, jph) + error = meteorsArray1[:,-1] + ind1 = numpy.where(error==0)[0] + penalty[ix,iy] = ind1.size + + i,j = numpy.unravel_index(penalty.argmax(), penalty.shape) + phOffset = jph_array[:,i,j] + + center_xangle = phOffset[pairx[1]] + center_yangle = phOffset[pairy[1]] + + phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j])) + phOffset = phOffset*180/numpy.pi + return phOffset + + + def run(self, dataOut, pairs, distances, hmin, hmax, nHours = None): + + dataOut.flagNoData = True + self.__dataReady = False + + if nHours == None: + nHours = 1 + + dataOut.outputInterval = nHours*3600 + + if self.__isConfig == False: +# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) + #Get Initial LTC time + self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) + self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() + + self.__isConfig = True + + if self.__buffer == None: + self.__buffer = dataOut.data_param.copy() + + else: + self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param)) + + self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready + + if self.__dataReady: + dataOut.utctimeInit = self.__initime + self.__initime += dataOut.outputInterval #to erase time offset + + freq = dataOut.frequency + c = dataOut.C #m/s + lamb = c/freq + k = 2*numpy.pi/lamb + azimuth = 0 + h = (hmin, hmax) + pairsList = ((0,3),(1,2)) + + meteorsArray = self.__buffer[0,:,:] + error = meteorsArray[:,-1] + boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14) + ind1 = numpy.where(boolError)[0] + meteorsArray = meteorsArray[ind1,:] + meteorsArray[:,-1] = 0 + phases = meteorsArray[:,9:13] + + #Calculate Gammas + gammas = self.__getGammas(pairs, k, distances, phases) +# gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180 + #Calculate Phases + phasesOff = self.__getPhases(azimuth, h, pairsList, distances, gammas, meteorsArray) + phasesOff = phasesOff.reshape((1,phasesOff.size)) + dataOut.data_output = -phasesOff + dataOut.flagNoData = False + self.__buffer = None + + + return + +class MeteorOperations(): + + def __init__(self): + + return + + def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, jph): + + arrayParameters = arrayParameters0.copy() + hmin = h[0] + hmax = h[1] + + #Calculate AOA (Error N 3, 4) + #JONES ET AL. 1998 + AOAthresh = numpy.pi/8 + error = arrayParameters[:,-1] + phases = -arrayParameters[:,9:13] + jph + arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth) + + #Calculate Heights (Error N 13 and 14) + error = arrayParameters[:,-1] + Ranges = arrayParameters[:,2] + zenith = arrayParameters[:,5] + arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax) + + #----------------------- Get Final data ------------------------------------ +# error = arrayParameters[:,-1] +# ind1 = numpy.where(error==0)[0] +# arrayParameters = arrayParameters[ind1,:] + + return arrayParameters + + def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth): + + arrayAOA = numpy.zeros((phases.shape[0],3)) + cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList) + + arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth) + cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1) + arrayAOA[:,2] = cosDirError + azimuthAngle = arrayAOA[:,0] + zenithAngle = arrayAOA[:,1] + + #Setting Error + #Number 3: AOA not fesible + indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0] + error[indInvalid] = 3 + #Number 4: Large difference in AOAs obtained from different antenna baselines + indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0] + error[indInvalid] = 4 + return arrayAOA, error + + def __getDirectionCosines(self, arrayPhase, pairsList): + #Initializing some variables + ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi + ang_aux = ang_aux.reshape(1,ang_aux.size) + + cosdir = numpy.zeros((arrayPhase.shape[0],2)) + cosdir0 = numpy.zeros((arrayPhase.shape[0],2)) + + + for i in range(2): + #First Estimation + phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]] + #Dealias + indcsi = numpy.where(phi0_aux > numpy.pi) + phi0_aux[indcsi] -= 2*numpy.pi + indcsi = numpy.where(phi0_aux < -numpy.pi) + phi0_aux[indcsi] += 2*numpy.pi + #Direction Cosine 0 + cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5) + + #Most-Accurate Second Estimation + phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]] + phi1_aux = phi1_aux.reshape(phi1_aux.size,1) + #Direction Cosine 1 + cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5) + + #Searching the correct Direction Cosine + cosdir0_aux = cosdir0[:,i] + cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1) + #Minimum Distance + cosDiff = (cosdir1 - cosdir0_aux)**2 + indcos = cosDiff.argmin(axis = 1) + #Saving Value obtained + cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos] + + return cosdir0, cosdir + + def __calculateAOA(self, cosdir, azimuth): + cosdirX = cosdir[:,0] + cosdirY = cosdir[:,1] + + zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi + azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east + angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose() + + return angles + + def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight): + + Ramb = 375 #Ramb = c/(2*PRF) + Re = 6371 #Earth Radius + heights = numpy.zeros(Ranges.shape) + + R_aux = numpy.array([0,1,2])*Ramb + R_aux = R_aux.reshape(1,R_aux.size) + + Ranges = Ranges.reshape(Ranges.size,1) + + Ri = Ranges + R_aux + hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re + + #Check if there is a height between 70 and 110 km + h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1) + ind_h = numpy.where(h_bool == 1)[0] + + hCorr = hi[ind_h, :] + ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight)) + + hCorr = hi[ind_hCorr] + heights[ind_h] = hCorr + + #Setting Error + #Number 13: Height unresolvable echo: not valid height within 70 to 110 km + #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km + + indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0] + error[indInvalid2] = 14 + indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0] + error[indInvalid1] = 13 + + return heights, error + \ No newline at end of file diff --git a/schainpy/scripts/Meteor_JASMET30MHz.py b/schainpy/scripts/Meteor_JASMET30MHz.py new file mode 100644 index 0000000..86d00ce --- /dev/null +++ b/schainpy/scripts/Meteor_JASMET30MHz.py @@ -0,0 +1,188 @@ +# DIAS 19 Y 20 FEB 2014 +# Comprobacion de Resultados DBS con SA + +import os, sys + +path = os.path.split(os.getcwd())[0] +path = os.path.split(path)[0] + +sys.path.insert(0, path) + +from schainpy.controller import Project + +desc = "JASMET Experiment Test" +filename = "JASMETtest.xml" + +controllerObj = Project() + +controllerObj.setup(id = '191', name='test01', description=desc) + +#Experimentos + +#2014051 20 Feb 2014 +# path = '/home/soporte/Data/JASMET/JASMET_30/2014106' +# pathFigure = '/home/soporte/workspace/Graficos/JASMET/prueba1' + +remotefolder = "/home/wmaster/graficos" +path = '/media/joscanoa/84A65E64A65E5730/soporte/Data/JASMET/JASMET_30/' + +pathfile1 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/meteor') +pathfile2 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/wind') +pathfile3 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/phase') + +pathfig = os.path.join(os.environ['HOME'],'Pictures/JASMET30/graph') + +startTime = '00:00:00' +endTime = '23:59:59' +xmin ='17.0' +xmax = '34.0' + +#------------------------------------------------------------------------------------------------ +readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader', + path=path, + startDate='2014/04/15', + endDate='2014/04/16', + startTime=startTime, + endTime=endTime, + online=0, + delay=5, + walk=1) + +opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') + + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId()) + +opObj00 = procUnitConfObj0.addOperation(name='selectChannels') +opObj00.addParameter(name='channelList', value='0, 1, 2, 3, 4', format='intlist') + +opObj01 = procUnitConfObj0.addOperation(name='setRadarFrequency') +# opObj01.addParameter(name='frequency', value='30.e6', format='float') +opObj01.addParameter(name='frequency', value='50.e6', format='float') + +#opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj1 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj0.getId()) +procUnitConfObj1.addParameter(name='nSeconds', value='100', format='int') + +opObj10 = procUnitConfObj1.addOperation(name='MeteorDetection') +# opObj10.addParameter(name='predefinedPhaseShifts', value='-89.5, 41.5, 0.0, -138.0, -85.5', format='floatlist') +opObj10.addParameter(name='predefinedPhaseShifts', value='0, 0, 0, 0, 0', format='floatlist') +opObj10.addParameter(name='cohDetection', value='0', format='bool') +opObj10.addParameter(name='noise_multiple', value='4', format='int') +opObj10.addParameter(name='SNRThresh', value='5', format='float') +opObj10.addParameter(name='phaseThresh', value='20', format='float') +opObj10.addParameter(name='azimuth', value='45', format='float') +opObj10.addParameter(name='hmin', value='68', format='float') +opObj10.addParameter(name='hmax', value='112', format='float') +opObj10.addParameter(name='saveAll', value='1', format='bool') + +opObj12 = procUnitConfObj1.addOperation(name='HDF5Writer', optype='other') +opObj12.addParameter(name='path', value=pathfile1) +opObj12.addParameter(name='blocksPerFile', value='1000', format='int') +opObj12.addParameter(name='metadataList',value='type,inputUnit,heightList,paramInterval',format='list') +opObj12.addParameter(name='dataList',value='data_param',format='list') +opObj12.addParameter(name='mode',value='0',format='int') +#Tiene que ser de 3 dimensiones, append en lugar de aumentar una dimension + +opObj13 = procUnitConfObj1.addOperation(name='SkyMapPlot', optype='other') +opObj13.addParameter(name='id', value='1', format='int') +opObj13.addParameter(name='wintitle', value='Sky Map', format='str') +opObj13.addParameter(name='save', value='1', format='bool') +opObj13.addParameter(name='figpath', value=pathfig, format='str') +opObj13.addParameter(name='ftp', value='1', format='int') +opObj13.addParameter(name='exp_code', value='15', format='int') +opObj13.addParameter(name='sub_exp_code', value='1', format='int') + + +#-------------------------------------------------------------------------------------------------- +procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId()) + +opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other') +opObj22.addParameter(name='technique', value='Meteors', format='str') +opObj22.addParameter(name='nHours', value='0.25', format='float') +opObj22.addParameter(name='SNRThresh', value='12.0', format='float') + +opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other') +opObj23.addParameter(name='id', value='2', format='int') +opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str') +opObj23.addParameter(name='save', value='1', format='bool') +opObj23.addParameter(name='figpath', value = pathfig, format='str') +opObj23.addParameter(name='zmin', value='-120', format='int') +opObj23.addParameter(name='zmax', value='120', format='int') +# opObj12.addParameter(name='zmin_ver', value='-0.8', format='float') +# opObj12.addParameter(name='zmax_ver', value='0.8', format='float') +# opObj23.addParameter(name='SNRmin', value='-10', format='int') +# opObj23.addParameter(name='SNRmax', value='60', format='int') +# opObj23.addParameter(name='SNRthresh', value='0', format='float') +opObj23.addParameter(name='xmin', value=xmin, format='float') +opObj23.addParameter(name='xmax', value=xmax, format='float') +opObj23.addParameter(name='ftp', value='1', format='int') +opObj23.addParameter(name='exp_code', value='15', format='int') +opObj23.addParameter(name='sub_exp_code', value='1', format='int') + +opObj24 = procUnitConfObj2.addOperation(name='HDF5Writer', optype='other') +opObj24.addParameter(name='path', value=pathfile2) +opObj24.addParameter(name='blocksPerFile', value='1000', format='int') +opObj24.addParameter(name='metadataList',value='type,inputUnit,outputInterval',format='list') +opObj24.addParameter(name='dataList',value='data_output,utctime',format='list') + +#-------------------------------------------------------------------------------------------------- +procUnitConfObj3 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId()) + + +opObj31 = procUnitConfObj3.addOperation(name='PhaseCalibration', optype='other') +opObj31.addParameter(name='nHours', value='0.25', format='float') +opObj31.addParameter(name='distances', value='-15, -15, 12, 12', format='intlist') +# opObj31.addParameter(name='distances', value='-25, -25, 20, 20', format='intlist') +opObj31.addParameter(name='pairs', value='(0,3),(1,2)', format='pairslist') +opObj31.addParameter(name='hmin', value='68', format='float') +opObj31.addParameter(name='hmax', value='112', format='float') + +opObj32 = procUnitConfObj3.addOperation(name='PhasePlot', optype='other') +opObj32.addParameter(name='id', value='201', format='int') +opObj32.addParameter(name='wintitle', value='PhaseCalibration', format='str') +# opObj32.addParameter(name='timerange', value='300', format='int') +opObj32.addParameter(name='xmin', value=xmin, format='float') +opObj32.addParameter(name='xmax', value=xmax, format='float') +# opObj32.addParameter(name='xmin', value='19', format='float') +# opObj32.addParameter(name='xmax', value='20', format='float') +opObj32.addParameter(name='ymin', value='-180', format='float') +opObj32.addParameter(name='ymax', value='180', format='float') +opObj32.addParameter(name='figpath', value=pathfig, format='str') +opObj32.addParameter(name='ftp', value='1', format='int') +opObj32.addParameter(name='exp_code', value='15', format='int') +opObj32.addParameter(name='sub_exp_code', value='1', format='int') + +opObj33 = procUnitConfObj3.addOperation(name='HDF5Writer', optype='other') +opObj33.addParameter(name='path', value=pathfile3) +opObj33.addParameter(name='blocksPerFile', value='1000', format='int') +opObj33.addParameter(name='metadataList',value='type,inputUnit,outputInterval',format='list') +opObj33.addParameter(name='dataList',value='data_output,utctime',format='list') +# opObj25.addParameter(name='mode',value='1,0,0',format='intlist') + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj4 = controllerObj.addProcUnit(name='SendToServer') +procUnitConfObj4.addParameter(name='server', value='jro-app.igp.gob.pe', format='str') +procUnitConfObj4.addParameter(name='username', value='wmaster', format='str') +procUnitConfObj4.addParameter(name='password', value='mst2010vhf', format='str') +procUnitConfObj4.addParameter(name='localfolder', value=pathfig, format='str') +procUnitConfObj4.addParameter(name='remotefolder', value=remotefolder, format='str') +procUnitConfObj4.addParameter(name='ext', value='.png', format='str') +procUnitConfObj4.addParameter(name='period', value=20, format='int') +procUnitConfObj4.addParameter(name='protocol', value='ftp', format='str') + +#-------------------------------------------------------------------------------------------------- + +print "Escribiendo el archivo XML" +controllerObj.writeXml(filename) +print "Leyendo el archivo XML" +controllerObj.readXml(filename) + +controllerObj.createObjects() +controllerObj.connectObjects() +controllerObj.run() \ No newline at end of file diff --git a/schainpy/scripts/Meteor_JASMET_PreProc.py b/schainpy/scripts/Meteor_JASMET_PreProc.py new file mode 100644 index 0000000..4b5d668 --- /dev/null +++ b/schainpy/scripts/Meteor_JASMET_PreProc.py @@ -0,0 +1,56 @@ +""" +Se debe verficar que el disco de datos se encuentra montado en el sistema +""" +import os, sys + +path = os.path.split(os.getcwd())[0] +path = os.path.split(path)[0] + +sys.path.insert(0, path) + +from schainpy.controller import Project + +desc = "Meteor Experiment Test" +filename = "meteor20130812.xml" + +controllerObj = Project() +controllerObj.setup(id = '191', name='meteor_test01', description=desc) + +path = '/home/dsuarez/.gvfs/data on 10.10.20.13/Jasmet50' + + +readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage', + path=path, + startDate='2014/04/16', + endDate='2014/04/16', + startTime='00:00:00', + endTime='23:59:59', + online=0, + walk=1) + +opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') + +procUnitConfObj0 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId()) + + +opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') + + +opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other') +opObj11.addParameter(name='n', value='2', format='int') + +opObj11 = procUnitConfObj0.addOperation(name='VoltageWriter', optype='other') +opObj11.addParameter(name='path', value='/home/jasmet/jasmet30_abril') +opObj11.addParameter(name='blocksPerFile', value='100', format='int') +opObj11.addParameter(name='profilesPerBlock', value='200', format='int') + + + +print "Escribiendo el archivo XML" +controllerObj.writeXml(filename) +print "Leyendo el archivo XML" +controllerObj.readXml(filename) + +controllerObj.createObjects() +controllerObj.connectObjects() +controllerObj.run()