From f590e95f1fbcf74c9837eee26f4c285cb00438d3 2020-04-21 21:03:35 From: avaldez Date: 2020-04-21 21:03:35 Subject: [PATCH] Update de WR-Project --- diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index dbaea8f..1040be4 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -548,7 +548,7 @@ class Spectra(JROData): deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor) velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - + if self.nmodes: return velrange/self.nmodes else: @@ -1104,7 +1104,7 @@ class PlotterData(object): MAXNUMY = 100 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False): - + self.key = code self.throttle = throttle_value self.exp_code = exp_code @@ -1139,7 +1139,7 @@ class PlotterData(object): return len(self.__times) def __getitem__(self, key): - + if key not in self.data: raise KeyError(log.error('Missing key: {}'.format(key))) if 'spc' in key or not self.buffering: @@ -1172,7 +1172,7 @@ class PlotterData(object): elif 'spc_moments' == plot: plot = 'moments' self.data[plot] = {} - + if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data: self.data['noise'] = {} self.data['rti'] = {} @@ -1180,7 +1180,7 @@ class PlotterData(object): self.plottypes.append('noise') if 'rti' not in self.plottypes: self.plottypes.append('rti') - + def shape(self, key): ''' Get the shape of the one-element data for the given key @@ -1196,17 +1196,17 @@ class PlotterData(object): ''' Update data object with new dataOut ''' - + if tm in self.__times: return self.profileIndex = dataOut.profileIndex self.tm = tm self.type = dataOut.type self.parameters = getattr(dataOut, 'parameters', []) - + if hasattr(dataOut, 'meta'): self.meta.update(dataOut.meta) - + self.pairs = dataOut.pairsList self.interval = dataOut.getTimeInterval() self.localtime = dataOut.useLocalTime @@ -1217,7 +1217,7 @@ class PlotterData(object): self.__heights.append(dataOut.heightList) self.__all_heights.update(dataOut.heightList) self.__times.append(tm) - + for plot in self.plottypes: if plot in ('spc', 'spc_moments'): z = dataOut.data_spc/dataOut.normFactor @@ -1250,8 +1250,8 @@ class PlotterData(object): if plot == 'scope': buffer = dataOut.data self.flagDataAsBlock = dataOut.flagDataAsBlock - self.nProfiles = dataOut.nProfiles - + self.nProfiles = dataOut.nProfiles + if plot == 'spc': self.data['spc'] = buffer elif plot == 'cspc': @@ -1326,7 +1326,7 @@ class PlotterData(object): else: meta['xrange'] = [] - meta.update(self.meta) + meta.update(self.meta) ret['metadata'] = meta return json.dumps(ret) diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index c3a39fe..da589d6 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -42,7 +42,7 @@ class SpectraPlot_(Figure): self.__xfilter_ena = False self.__yfilter_ena = False - + self.indice=1 def getSubplots(self): @@ -225,11 +225,231 @@ class SpectraPlot_(Figure): ftp=ftp, wr_period=wr_period, thisDatetime=thisDatetime) - + return dataOut @MPDecorator +class WpowerPlot_(Figure): + + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'wpo' + + def __init__(self): + Figure.__init__(self) + 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 = WPO_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 + + self.indice=1 + + 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): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : None + """ + print("***************PLOTEO******************") + print("DATAOUT SHAPE : ",dataOut.data.shape) + if dataOut.flagNoData: + return dataOut + + 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)) + + + print("channelIndexList",channelIndexList) + if normFactor is None: + factor = dataOut.normFactor + else: + factor = normFactor + if xaxis == "frequency": + x = dataOut.getFreqRange(1)/1000. + xlabel = "Frequency (kHz)" + + elif xaxis == "time": + x = dataOut.getAcfRange(1) + xlabel = "Time (ms)" + + else: + x = dataOut.getVelRange(1) + xlabel = "Velocity (m/s)" + + ylabel = "Range (km)" + + y = dataOut.getHeiRange() + print("factor",factor) + + z = dataOut.data/factor # dividido /factor + 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.getNoise()/factor + noisedB = 10*numpy.log10(noise) + + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) + title = wintitle + "Weather Power" + + 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) + return dataOut + +@MPDecorator class CrossSpectraPlot_(Figure): isConfig = None @@ -256,7 +476,7 @@ class CrossSpectraPlot_(Figure): self.EXP_CODE = None self.SUB_EXP_CODE = None self.PLOT_POS = None - + self.indice=0 def getSubplots(self): @@ -314,7 +534,7 @@ class CrossSpectraPlot_(Figure): zmax : None """ - if dataOut.flagNoData: + if dataOut.flagNoData: return dataOut if pairsList == None: @@ -331,7 +551,7 @@ class CrossSpectraPlot_(Figure): if len(pairsIndexList) > 4: pairsIndexList = pairsIndexList[0:4] - + if normFactor is None: factor = dataOut.normFactor else: @@ -402,7 +622,7 @@ class CrossSpectraPlot_(Figure): self.isConfig = True self.setWinTitle(title) - + for i in range(self.nplots): pair = dataOut.pairsList[pairsIndexList[i]] @@ -563,7 +783,7 @@ class RTIPlot_(Figure): #colormap = kwargs.get('colormap', 'jet') if HEIGHT is not None: self.HEIGHT = HEIGHT - + if not isTimeInHourRange(dataOut.datatime, xmin, xmax): return @@ -745,7 +965,7 @@ class CoherenceMap_(Figure): ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): - if dataOut.flagNoData: + if dataOut.flagNoData: return dataOut if not isTimeInHourRange(dataOut.datatime, xmin, xmax): @@ -935,7 +1155,7 @@ class PowerProfilePlot_(Figure): ftp=False, wr_period=1, server=None, folder=None, username=None, password=None): - if dataOut.flagNoData: + if dataOut.flagNoData: return dataOut @@ -1009,7 +1229,7 @@ class PowerProfilePlot_(Figure): ftp=ftp, wr_period=wr_period, thisDatetime=thisDatetime) - + return dataOut @MPDecorator @@ -1066,7 +1286,7 @@ class SpectraCutPlot_(Figure): folder=None, username=None, password=None, xaxis="frequency"): - if dataOut.flagNoData: + if dataOut.flagNoData: return dataOut if channelList == None: @@ -1248,7 +1468,7 @@ class Noise_(Figure): server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): - if dataOut.flagNoData: + if dataOut.flagNoData: return dataOut if not isTimeInHourRange(dataOut.datatime, xmin, xmax): @@ -1444,7 +1664,7 @@ class BeaconPhase_(Figure): server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): - if dataOut.flagNoData: + if dataOut.flagNoData: return dataOut if not isTimeInHourRange(dataOut.datatime, xmin, xmax): @@ -1586,4 +1806,4 @@ class BeaconPhase_(Figure): thisDatetime=thisDatetime, update_figfile=update_figfile) - return dataOut \ No newline at end of file + return dataOut diff --git a/schainpy/model/graphics/jroplot_voltage.py b/schainpy/model/graphics/jroplot_voltage.py index ae4a972..be5c0c1 100644 --- a/schainpy/model/graphics/jroplot_voltage.py +++ b/schainpy/model/graphics/jroplot_voltage.py @@ -13,62 +13,62 @@ from .figure import Figure @MPDecorator class Scope_(Figure): - + isConfig = None - + def __init__(self):#, **kwargs): #YONG Figure.__init__(self)#, **kwargs) self.isConfig = False self.WIDTH = 300 self.HEIGHT = 200 self.counter_imagwr = 0 - + def getSubplots(self): - + nrow = self.nplots ncol = 3 return nrow, ncol - + def setup(self, id, nplots, wintitle, show): - + self.nplots = nplots - - self.createFigure(id=id, - wintitle=wintitle, + + self.createFigure(id=id, + wintitle=wintitle, show=show) - + nrow,ncol = self.getSubplots() colspan = 3 rowspan = 1 - + for i in range(nplots): self.addAxes(nrow, ncol, i, 0, colspan, rowspan) - + def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): yreal = y[channelIndexList,:].real yimag = y[channelIndexList,:].imag - + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) xlabel = "Range (Km)" ylabel = "Intensity - IQ" - + if not self.isConfig: nplots = len(channelIndexList) - + self.setup(id=id, nplots=nplots, wintitle='', show=show) - + if xmin == None: xmin = numpy.nanmin(x) if xmax == None: xmax = numpy.nanmax(x) if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag)) if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag)) - + self.isConfig = True - + self.setWinTitle(title) - + for i in range(len(self.axesList)): title = "Channel %d" %(i) axes = self.axesList[i] @@ -78,32 +78,32 @@ class Scope_(Figure): xlabel=xlabel, ylabel=ylabel, title=title) axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2) - + def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:]) yreal = y.real - + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) xlabel = "Range (Km)" ylabel = "Intensity" - + if not self.isConfig: nplots = len(channelIndexList) - + self.setup(id=id, nplots=nplots, wintitle='', show=show) - + if xmin == None: xmin = numpy.nanmin(x) if xmax == None: xmax = numpy.nanmax(x) if ymin == None: ymin = numpy.nanmin(yreal) if ymax == None: ymax = numpy.nanmax(yreal) - + self.isConfig = True - + self.setWinTitle(title) - + for i in range(len(self.axesList)): title = "Channel %d" %(i) axes = self.axesList[i] @@ -112,14 +112,48 @@ class Scope_(Figure): xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, xlabel=xlabel, ylabel=ylabel, title=title) - + def plot_weatherpower(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + y = y[channelIndexList,:] + yreal = y + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(yreal) + if ymax == None: ymax = numpy.nanmax(yreal) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + ychannel = yreal[i,:] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + + def run(self, dataOut, id, wintitle="", channelList=None, xmin=None, xmax=None, ymin=None, ymax=None, save=False, figpath='./', figfile=None, show=True, wr_period=1, ftp=False, server=None, folder=None, username=None, password=None, type='power', **kwargs): - + """ - + Input: dataOut : id : @@ -130,9 +164,9 @@ class Scope_(Figure): ymin : None, ymax : None, """ - if dataOut.flagNoData: + if dataOut.flagNoData: return dataOut - + if channelList == None: channelIndexList = dataOut.channelIndexList else: @@ -141,20 +175,22 @@ class Scope_(Figure): if channel not in dataOut.channelList: raise ValueError("Channel %d is not in dataOut.channelList") channelIndexList.append(dataOut.channelList.index(channel)) - + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) - + ### print("***************** PLOTEO **************************") + ### print(dataOut.nProfiles) + ### print(dataOut.heightList.shape) if dataOut.flagDataAsBlock: - + for i in range(dataOut.nProfiles): - + wintitle1 = wintitle + " [Profile = %d] " %i - + if type == "power": - self.plot_power(dataOut.heightList, + self.plot_power(dataOut.heightList, dataOut.data[:,i,:], - id, - channelIndexList, + id, + channelIndexList, thisDatetime, wintitle1, show, @@ -162,12 +198,38 @@ class Scope_(Figure): xmax, ymin, ymax) - + + if type == "weatherpower": + self.plot_weatherpower(dataOut.heightList, + dataOut.data[:,i,:], + id, + channelIndexList, + thisDatetime, + wintitle1, + show, + xmin, + xmax, + ymin, + ymax) + + if type == "weathervelocity": + self.plot_weatherpower(dataOut.heightList, + dataOut.data_velocity[:,i,:], + id, + channelIndexList, + thisDatetime, + wintitle1, + show, + xmin, + xmax, + ymin, + ymax) + if type == "iq": - self.plot_iq(dataOut.heightList, + self.plot_iq(dataOut.heightList, dataOut.data[:,i,:], - id, - channelIndexList, + id, + channelIndexList, thisDatetime, wintitle1, show, @@ -175,27 +237,27 @@ class Scope_(Figure): xmax, ymin, ymax) - + self.draw() - + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") figfile = self.getFilename(name = str_datetime) + "_" + str(i) - + self.save(figpath=figpath, figfile=figfile, save=save, ftp=ftp, wr_period=wr_period, thisDatetime=thisDatetime) - + else: wintitle += " [Profile = %d] " %dataOut.profileIndex - + if type == "power": - self.plot_power(dataOut.heightList, + self.plot_power(dataOut.heightList, dataOut.data, - id, - channelIndexList, + id, + channelIndexList, thisDatetime, wintitle, show, @@ -203,12 +265,12 @@ class Scope_(Figure): xmax, ymin, ymax) - + if type == "iq": - self.plot_iq(dataOut.heightList, + self.plot_iq(dataOut.heightList, dataOut.data, - id, - channelIndexList, + id, + channelIndexList, thisDatetime, wintitle, show, @@ -216,12 +278,12 @@ class Scope_(Figure): xmax, ymin, ymax) - + self.draw() - + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + "_" + str(dataOut.profileIndex) - figfile = self.getFilename(name = str_datetime) - + figfile = self.getFilename(name = str_datetime) + self.save(figpath=figpath, figfile=figfile, save=save, @@ -229,4 +291,4 @@ class Scope_(Figure): wr_period=wr_period, thisDatetime=thisDatetime) - return dataOut \ No newline at end of file + return dataOut diff --git a/schainpy/model/graphics/plotting_codes.py b/schainpy/model/graphics/plotting_codes.py index 3432f0c..52bd98d 100644 --- a/schainpy/model/graphics/plotting_codes.py +++ b/schainpy/model/graphics/plotting_codes.py @@ -23,6 +23,8 @@ MSKYMAP_CODE = 23 MPHASE_CODE = 24 MOMENTS_CODE = 25 -PARMS_CODE = 26 +PARMS_CODE = 26 SPECFIT_CODE = 27 EWDRIFT_CODE = 28 + +WPO_CODE = 29 #Weather Intensity - Power diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index cdb91da..2d91134 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -183,7 +183,7 @@ class ParamReader(JRODataReader,ProcessingUnit): except IOError: traceback.print_exc() raise IOError("The file %s can't be opened" %(filename)) - + #In case has utctime attribute grp2 = grp1['utctime'] # thisUtcTime = grp2.value[0] - 5*3600 #To convert to local time @@ -497,7 +497,7 @@ class ParamWriter(Operation): setType = None def __init__(self): - + Operation.__init__(self) return @@ -530,9 +530,9 @@ class ParamWriter(Operation): dsDict['variable'] = self.dataList[i] #--------------------- Conditionals ------------------------ #There is no data - + if dataAux is None: - + return 0 if isinstance(dataAux, (int, float, numpy.integer, numpy.float)): @@ -704,7 +704,7 @@ class ParamWriter(Operation): return False def setNextFile(self): - + ext = self.ext path = self.path setFile = self.setFile @@ -785,7 +785,7 @@ class ParamWriter(Operation): for j in range(dsInfo['dsNumber']): dsInfo = dsList[i] tableName = dsInfo['dsName'] - + if dsInfo['nDim'] == 3: shape = dsInfo['shape'].astype(int) @@ -869,7 +869,8 @@ class ParamWriter(Operation): dsList = self.dsList data = self.data ind = 0 - + #print("dsList ",dsList) + #print("len ",len(dsList)) while ind < len(dsList): dsInfo = dsList[ind] dataAux = getattr(self.dataOut, dsInfo['variable']) @@ -902,22 +903,41 @@ class ParamWriter(Operation): dsList = self.dsList for i in range(len(self.ds)): + print("#############", i , "#######################") dsInfo = dsList[i] nDim = dsInfo['nDim'] mode = dsInfo['mode'] - + print("dsInfo",dsInfo) + print("nDim",nDim) + print("mode",mode) # First time if self.firsttime: + print("ENTRE FIRSTIME") if type(self.data[i]) == numpy.ndarray: if nDim == 3: + print("xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx") + print("ndim","dentro del primer if 3") self.data[i] = self.data[i].reshape((self.data[i].shape[0],self.data[i].shape[1],1)) + print(self.data[i].shape) + print(type(self.data[i])) self.ds[i].resize(self.data[i].shape) + print(self.ds[i].shape) + print(type(self.ds[i])) if mode == 2: self.ds[i].resize(self.data[i].shape) - self.ds[i][:] = self.data[i] - else: + try: + print("PTM ODIO ESTO") + print(self.ds[i][:].shape) + self.ds[i][:] = self.data[i] + print("*****___________********______******") + except: + print("q habra pasaado") + return + print("LLEGUE Y CUMPLI EL IF") + else: + print("ELSE -----------------------") # From second time # Meteors! if mode == 2: @@ -944,7 +964,8 @@ class ParamWriter(Operation): self.firsttime = False self.blockIndex += 1 - + print("HOLA AMIGOS COMO ESTAN LLEGUE") + print("xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx") #Close to save changes self.fp.flush() self.fp.close() @@ -954,7 +975,7 @@ class ParamWriter(Operation): self.dataOut = dataOut if not(self.isConfig): - self.setup(dataOut, path=path, blocksPerFile=blocksPerFile, + self.setup(dataOut, path=path, blocksPerFile=blocksPerFile, metadataList=metadataList, dataList=dataList, mode=mode, setType=setType) @@ -963,7 +984,7 @@ class ParamWriter(Operation): self.putData() return - + @MPDecorator class ParameterReader(Reader, ProcessingUnit): @@ -992,43 +1013,43 @@ class ParameterReader(Reader, ProcessingUnit): self.set_kwargs(**kwargs) if not self.ext.startswith('.'): - self.ext = '.{}'.format(self.ext) + self.ext = '.{}'.format(self.ext) if self.online: log.log("Searching files in online mode...", self.name) for nTries in range(self.nTries): fullpath = self.searchFilesOnLine(self.path, self.startDate, - self.endDate, self.expLabel, self.ext, self.walk, + self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt) try: fullpath = next(fullpath) except: fullpath = None - + if fullpath: break log.warning( 'Waiting {} sec for a valid file in {}: try {} ...'.format( - self.delay, self.path, nTries + 1), + self.delay, self.path, nTries + 1), self.name) time.sleep(self.delay) if not(fullpath): raise schainpy.admin.SchainError( - 'There isn\'t any valid file in {}'.format(self.path)) + 'There isn\'t any valid file in {}'.format(self.path)) pathname, filename = os.path.split(fullpath) self.year = int(filename[1:5]) self.doy = int(filename[5:8]) - self.set = int(filename[8:11]) - 1 + self.set = int(filename[8:11]) - 1 else: log.log("Searching files in {}".format(self.path), self.name) - self.filenameList = self.searchFilesOffLine(self.path, self.startDate, + self.filenameList = self.searchFilesOffLine(self.path, self.startDate, self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt) - + self.setNextFile() return @@ -1036,11 +1057,11 @@ class ParameterReader(Reader, ProcessingUnit): def readFirstHeader(self): '''Read metadata and data''' - self.__readMetadata() + self.__readMetadata() self.__readData() self.__setBlockList() self.blockIndex = 0 - + return def __setBlockList(self): @@ -1099,7 +1120,7 @@ class ParameterReader(Reader, ProcessingUnit): else: data = gp[name].value listMetaname.append(name) - listMetadata.append(data) + listMetadata.append(data) elif self.metadata: metadata = json.loads(self.metadata) listShapes = {} @@ -1115,7 +1136,7 @@ class ParameterReader(Reader, ProcessingUnit): self.listShapes = listShapes self.listMetaname = listMetaname - self.listMeta = listMetadata + self.listMeta = listMetadata return @@ -1123,7 +1144,7 @@ class ParameterReader(Reader, ProcessingUnit): listdataname = [] listdata = [] - + if 'Data' in self.fp: grp = self.fp['Data'] for item in list(grp.items()): @@ -1137,7 +1158,7 @@ class ParameterReader(Reader, ProcessingUnit): for i in range(dim): array.append(grp[name]['table{:02d}'.format(i)].value) array = numpy.array(array) - + listdata.append(array) elif self.metadata: metadata = json.loads(self.metadata) @@ -1160,7 +1181,7 @@ class ParameterReader(Reader, ProcessingUnit): self.listDataname = listdataname self.listData = listdata return - + def getData(self): for i in range(len(self.listMeta)): @@ -1230,7 +1251,7 @@ class ParameterWriter(Operation): lastTime = None def __init__(self): - + Operation.__init__(self) return @@ -1257,7 +1278,7 @@ class ParameterWriter(Operation): dsDict['nDim'] = len(dataAux.shape) dsDict['shape'] = dataAux.shape dsDict['dsNumber'] = dataAux.shape[0] - + dsList.append(dsDict) tableList.append((self.dataList[i], dsDict['nDim'])) @@ -1274,7 +1295,7 @@ class ParameterWriter(Operation): self.lastTime = currentTime self.currentDay = dataDay return False - + timeDiff = currentTime - self.lastTime #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora @@ -1292,7 +1313,7 @@ class ParameterWriter(Operation): self.dataOut = dataOut if not(self.isConfig): - self.setup(path=path, blocksPerFile=blocksPerFile, + self.setup(path=path, blocksPerFile=blocksPerFile, metadataList=metadataList, dataList=dataList, setType=setType) @@ -1301,9 +1322,9 @@ class ParameterWriter(Operation): self.putData() return - + def setNextFile(self): - + ext = self.ext path = self.path setFile = self.setFile @@ -1369,17 +1390,17 @@ class ParameterWriter(Operation): return def writeData(self, fp): - + grp = fp.create_group("Data") dtsets = [] data = [] - + for dsInfo in self.dsList: if dsInfo['nDim'] == 0: ds = grp.create_dataset( - dsInfo['variable'], + dsInfo['variable'], (self.blocksPerFile, ), - chunks=True, + chunks=True, dtype=numpy.float64) dtsets.append(ds) data.append((dsInfo['variable'], -1)) @@ -1387,7 +1408,7 @@ class ParameterWriter(Operation): sgrp = grp.create_group(dsInfo['variable']) for i in range(dsInfo['dsNumber']): ds = sgrp.create_dataset( - 'table{:02d}'.format(i), + 'table{:02d}'.format(i), (self.blocksPerFile, ) + dsInfo['shape'][1:], chunks=True) dtsets.append(ds) @@ -1395,7 +1416,7 @@ class ParameterWriter(Operation): fp.flush() log.log('Creating file: {}'.format(fp.filename), self.name) - + self.ds = dtsets self.data = data self.firsttime = True @@ -1409,10 +1430,12 @@ class ParameterWriter(Operation): self.setNextFile() for i, ds in enumerate(self.ds): + print(i,ds) attr, ch = self.data[i] if ch == -1: ds[self.blockIndex] = getattr(self.dataOut, attr) else: + print(ch, getattr(self.dataOut, attr).shape) ds[self.blockIndex] = getattr(self.dataOut, attr)[ch] self.fp.flush() diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 571dfe7..d25a776 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -8,12 +8,12 @@ import copy import sys import importlib import itertools -from multiprocessing import Pool, TimeoutError +from multiprocessing import Pool, TimeoutError from multiprocessing.pool import ThreadPool import time from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters -from .jroproc_base import ProcessingUnit, Operation, MPDecorator +from .jroproc_base import ProcessingUnit, Operation, MPDecorator from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon from scipy import asarray as ar,exp from scipy.optimize import curve_fit @@ -48,13 +48,13 @@ def _unpickle_method(func_name, obj, cls): @MPDecorator class ParametersProc(ProcessingUnit): - + METHODS = {} nSeconds = None def __init__(self): ProcessingUnit.__init__(self) - + # self.objectDict = {} self.buffer = None self.firstdatatime = None @@ -63,14 +63,14 @@ class ParametersProc(ProcessingUnit): self.setupReq = False #Agregar a todas las unidades de proc def __updateObjFromInput(self): - + self.dataOut.inputUnit = self.dataIn.type - + self.dataOut.timeZone = self.dataIn.timeZone self.dataOut.dstFlag = self.dataIn.dstFlag self.dataOut.errorCount = self.dataIn.errorCount self.dataOut.useLocalTime = self.dataIn.useLocalTime - + self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy() self.dataOut.channelList = self.dataIn.channelList @@ -92,27 +92,32 @@ class ParametersProc(ProcessingUnit): self.dataOut.ippSeconds = self.dataIn.ippSeconds # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter self.dataOut.timeInterval1 = self.dataIn.timeInterval - self.dataOut.heightList = self.dataIn.getHeiRange() + self.dataOut.heightList = self.dataIn.getHeiRange() self.dataOut.frequency = self.dataIn.frequency # self.dataOut.noise = self.dataIn.noise - - def run(self): - - + def run(self): #---------------------- Voltage Data --------------------------- - if self.dataIn.type == "Voltage": - + print(" *************INSIDE PARAMETER********") self.__updateObjFromInput() self.dataOut.data_pre = self.dataIn.data.copy() self.dataOut.flagNoData = False self.dataOut.utctimeInit = self.dataIn.utctime - self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds + self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds + + if self.dataIn.identifierWR== True: + print("**********************************************") + self.dataOut.data_intensity = self.dataIn.data #valor para intensidad + self.dataOut.data_velocity = self.dataIn.data_velocity #valor para velocidad + self.dataOut.identifierWR = self.dataIn.identifierWR + self.dataOut.PRFbyAngle = self.dataIn.PRFbyAngle + print(self.dataOut.data_intensity.shape) + print(self.dataOut.utctimeInit) + print(self.dataOut.nCohInt ) + print(self.dataOut.PRFbyAngle) return - #---------------------- Spectra Data --------------------------- - if self.dataIn.type == "Spectra": self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc) @@ -126,243 +131,256 @@ class ParametersProc(ProcessingUnit): self.dataOut.spc_noise = self.dataIn.getNoise() self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1)) # self.dataOut.normFactor = self.dataIn.normFactor - self.dataOut.pairsList = self.dataIn.pairsList + self.dataOut.pairsList = self.dataIn.pairsList self.dataOut.groupList = self.dataIn.pairsList - self.dataOut.flagNoData = False - + self.dataOut.flagNoData = False + + if self.dataIn.flagWR== True: + print("##############################################") + self.dataOut.nIncohInt = self.dataIn.nIncohInt + self.dataOut.utctimeInit = self.dataIn.utctime + self.dataOut.data_intensity = self.dataIn.data #valor para intensidad + self.dataOut.data_velocity = self.dataIn.data_velocity #valor para velocidad + self.dataOut.flagWR = self.dataIn.flagWR + self.dataOut.PRFbyAngle = self.dataIn.PRFbyAngle + print(self.dataOut.data_intensity.shape) + print(self.dataOut.utctimeInit) + print(self.dataOut.nIncohInt ) + print(self.dataOut.PRFbyAngle) + if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels self.dataOut.ChanDist = self.dataIn.ChanDist - else: self.dataOut.ChanDist = None - + else: self.dataOut.ChanDist = None + #if hasattr(self.dataIn, 'VelRange'): #Velocities range # self.dataOut.VelRange = self.dataIn.VelRange #else: self.dataOut.VelRange = None - + if hasattr(self.dataIn, 'RadarConst'): #Radar Constant self.dataOut.RadarConst = self.dataIn.RadarConst - + if hasattr(self.dataIn, 'NPW'): #NPW self.dataOut.NPW = self.dataIn.NPW - + if hasattr(self.dataIn, 'COFA'): #COFA self.dataOut.COFA = self.dataIn.COFA - - - + + + #---------------------- Correlation Data --------------------------- - + if self.dataIn.type == "Correlation": acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions() - + self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:]) self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:]) self.dataOut.groupList = (acf_pairs, ccf_pairs) - + self.dataOut.abscissaList = self.dataIn.lagRange self.dataOut.noise = self.dataIn.noise self.dataOut.data_SNR = self.dataIn.SNR self.dataOut.flagNoData = False self.dataOut.nAvg = self.dataIn.nAvg - + #---------------------- Parameters Data --------------------------- - + if self.dataIn.type == "Parameters": self.dataOut.copy(self.dataIn) self.dataOut.flagNoData = False - + return True - + self.__updateObjFromInput() self.dataOut.utctimeInit = self.dataIn.utctime self.dataOut.paramInterval = self.dataIn.timeInterval - + return def target(tups): - + obj, args = tups - + return obj.FitGau(args) - - + + class SpectralFilters(Operation): - + '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR - + LimitR : It is the limit in m/s of Rainfall LimitW : It is the limit in m/s for Winds - + Input: - + self.dataOut.data_pre : SPC and CSPC self.dataOut.spc_range : To select wind and rainfall velocities - + Affected: - + self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind - self.dataOut.spcparam_range : Used in SpcParamPlot + self.dataOut.spcparam_range : Used in SpcParamPlot self.dataOut.SPCparam : Used in PrecipitationProc - - + + ''' - + def __init__(self): Operation.__init__(self) self.i=0 - - def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5): - - - #Limite de vientos + + def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5): + + + #Limite de vientos LimitR = PositiveLimit LimitN = NegativeLimit - + self.spc = dataOut.data_pre[0].copy() self.cspc = dataOut.data_pre[1].copy() - + self.Num_Hei = self.spc.shape[2] self.Num_Bin = self.spc.shape[1] self.Num_Chn = self.spc.shape[0] - + VelRange = dataOut.spc_range[2] TimeRange = dataOut.spc_range[1] FrecRange = dataOut.spc_range[0] - + Vmax= 2*numpy.max(dataOut.spc_range[2]) Tmax= 2*numpy.max(dataOut.spc_range[1]) Fmax= 2*numpy.max(dataOut.spc_range[0]) - + Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()] Breaker1R=numpy.where(VelRange == Breaker1R) - - Delta = self.Num_Bin/2 - Breaker1R[0] - - + + Delta = self.Num_Bin/2 - Breaker1R[0] + + '''Reacomodando SPCrange''' VelRange=numpy.roll(VelRange,-(int(self.Num_Bin/2)) ,axis=0) - + VelRange[-(int(self.Num_Bin/2)):]+= Vmax - + FrecRange=numpy.roll(FrecRange,-(int(self.Num_Bin/2)),axis=0) - + FrecRange[-(int(self.Num_Bin/2)):]+= Fmax - + TimeRange=numpy.roll(TimeRange,-(int(self.Num_Bin/2)),axis=0) - + TimeRange[-(int(self.Num_Bin/2)):]+= Tmax - + ''' ------------------ ''' - + Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()] Breaker2R=numpy.where(VelRange == Breaker2R) - - + + SPCroll = numpy.roll(self.spc,-(int(self.Num_Bin/2)) ,axis=1) - + SPCcut = SPCroll.copy() for i in range(self.Num_Chn): - + SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i] SPCcut[i,-int(Delta):,:] = dataOut.noise[i] - + SPCcut[i]=SPCcut[i]- dataOut.noise[i] SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20 - + SPCroll[i]=SPCroll[i]-dataOut.noise[i] SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20 - + SPC_ch1 = SPCroll - + SPC_ch2 = SPCcut - + SPCparam = (SPC_ch1, SPC_ch2, self.spc) - dataOut.SPCparam = numpy.asarray(SPCparam) - - + dataOut.SPCparam = numpy.asarray(SPCparam) + + dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1]) - + dataOut.spcparam_range[2]=VelRange dataOut.spcparam_range[1]=TimeRange dataOut.spcparam_range[0]=FrecRange return dataOut - + class GaussianFit(Operation): - + ''' - Function that fit of one and two generalized gaussians (gg) based - on the PSD shape across an "power band" identified from a cumsum of + Function that fit of one and two generalized gaussians (gg) based + on the PSD shape across an "power band" identified from a cumsum of the measured spectrum - noise. - + Input: self.dataOut.data_pre : SelfSpectra - + Output: self.dataOut.SPCparam : SPC_ch1, SPC_ch2 - + ''' def __init__(self): Operation.__init__(self) self.i=0 - - + + def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points """This routine will find a couple of generalized Gaussians to a power spectrum input: spc output: Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise """ - + self.spc = dataOut.data_pre[0].copy() self.Num_Hei = self.spc.shape[2] self.Num_Bin = self.spc.shape[1] self.Num_Chn = self.spc.shape[0] Vrange = dataOut.abscissaList - + GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei]) SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch1[:] = numpy.NaN SPC_ch2[:] = numpy.NaN - + start_time = time.time() - + noise_ = dataOut.spc_noise[0].copy() - - - pool = Pool(processes=self.Num_Chn) + + + pool = Pool(processes=self.Num_Chn) args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)] - objs = [self for __ in range(self.Num_Chn)] - attrs = list(zip(objs, args)) + objs = [self for __ in range(self.Num_Chn)] + attrs = list(zip(objs, args)) gauSPC = pool.map(target, attrs) dataOut.SPCparam = numpy.asarray(SPCparam) - + ''' Parameters: 1. Amplitude 2. Shift 3. Width 4. Power ''' - + def FitGau(self, X): - + Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X - + SPCparam = [] SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch1[:] = 0#numpy.NaN SPC_ch2[:] = 0#numpy.NaN - - - + + + for ht in range(self.Num_Hei): - - + + spc = numpy.asarray(self.spc)[ch,:,ht] - + ############################################# # normalizing spc and noise # This part differs from gg1 @@ -370,60 +388,60 @@ class GaussianFit(Operation): #spc = spc / spc_norm_max pnoise = pnoise #/ spc_norm_max ############################################# - + fatspectra=1.0 - + wnoise = noise_ #/ spc_norm_max #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used - #if wnoise>1.1*pnoise: # to be tested later + #if wnoise>1.1*pnoise: # to be tested later # wnoise=pnoise - noisebl=wnoise*0.9; + noisebl=wnoise*0.9; noisebh=wnoise*1.1 spc=spc-wnoise - + minx=numpy.argmin(spc) - #spcs=spc.copy() + #spcs=spc.copy() spcs=numpy.roll(spc,-minx) cum=numpy.cumsum(spcs) tot_noise=wnoise * self.Num_Bin #64; - + snr = sum(spcs)/tot_noise snrdB=10.*numpy.log10(snr) - + if snrdB < SNRlimit : snr = numpy.NaN SPC_ch1[:,ht] = 0#numpy.NaN SPC_ch1[:,ht] = 0#numpy.NaN SPCparam = (SPC_ch1,SPC_ch2) continue - - + + #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None - - cummax=max(cum); + + cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region - cumlo=cummax*epsi; + cumlo=cummax*epsi; cumhi=cummax*(1-epsi) powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum-12: # when SNR is strong pick the peak with least shift (LOS velocity) error if oneG: choice=0 @@ -481,10 +499,10 @@ class GaussianFit(Operation): w1=lsq2[0][1]; w2=lsq2[0][5] a1=lsq2[0][2]; a2=lsq2[0][6] p1=lsq2[0][3]; p2=lsq2[0][7] - s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; + s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2; gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling - + if gp1>gp2: if a1>0.7*a2: choice=1 @@ -499,157 +517,157 @@ class GaussianFit(Operation): choice=numpy.argmax([a1,a2])+1 #else: #choice=argmin([std2a,std2b])+1 - + else: # with low SNR go to the most energetic peak choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]]) - - - shift0=lsq2[0][0]; + + + shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0]) - shift1=lsq2[0][4]; + shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0]) - + max_vel = 1.0 - + #first peak will be 0, second peak will be 1 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range shift0=lsq2[0][0] width0=lsq2[0][1] Amplitude0=lsq2[0][2] p0=lsq2[0][3] - + shift1=lsq2[0][4] width1=lsq2[0][5] Amplitude1=lsq2[0][6] p1=lsq2[0][7] - noise=lsq2[0][8] + noise=lsq2[0][8] else: shift1=lsq2[0][0] width1=lsq2[0][1] Amplitude1=lsq2[0][2] p1=lsq2[0][3] - + shift0=lsq2[0][4] width0=lsq2[0][5] Amplitude0=lsq2[0][6] - p0=lsq2[0][7] - noise=lsq2[0][8] - + p0=lsq2[0][7] + noise=lsq2[0][8] + if Amplitude0<0.05: # in case the peak is noise - shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN] + shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN] if Amplitude1<0.05: - shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN] - - + shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN] + + SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1 SPCparam = (SPC_ch1,SPC_ch2) - - + + return GauSPC - + def y_model1(self,x,state): shift0,width0,amplitude0,power0,noise=state model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0) - + model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0) - + model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0) return model0+model0u+model0d+noise - - def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist + + def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0) - + model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0) - + model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0) model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1) - + model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1) - + model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1) return model0+model0u+model0d+model1+model1u+model1d+noise - - def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is. + + def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is. return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))**2)#/(64-5.) # /(64-5.) can be commented - + def misfit2(self,state,y_data,x,num_intg): return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.) - - + + class PrecipitationProc(Operation): - + ''' Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R) - - Input: + + Input: self.dataOut.data_pre : SelfSpectra - - Output: - - self.dataOut.data_output : Reflectivity factor, rainfall Rate - - - Parameters affected: + + Output: + + self.dataOut.data_output : Reflectivity factor, rainfall Rate + + + Parameters affected: ''' - + def __init__(self): Operation.__init__(self) self.i=0 - - + + def gaus(self,xSamples,Amp,Mu,Sigma): return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) )) - - - + + + def Moments(self, ySamples, xSamples): Pot = numpy.nansum( ySamples ) # Potencia, momento 0 yNorm = ySamples / Pot - + Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento - Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento + Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral - - return numpy.array([Pot, Vr, Desv]) - - def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118, + + return numpy.array([Pot, Vr, Desv]) + + def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118, tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350): - - + + Velrange = dataOut.spcparam_range[2] FrecRange = dataOut.spcparam_range[0] - + dV= Velrange[1]-Velrange[0] dF= FrecRange[1]-FrecRange[0] - + if radar == "MIRA35C" : - + self.spc = dataOut.data_pre[0].copy() self.Num_Hei = self.spc.shape[2] self.Num_Bin = self.spc.shape[1] self.Num_Chn = self.spc.shape[0] Ze = self.dBZeMODE2(dataOut) - + else: - + self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() # - + """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX""" - - self.spc[:,:,0:7]= numpy.NaN - + + self.spc[:,:,0:7]= numpy.NaN + """##########################################""" - + self.Num_Hei = self.spc.shape[2] self.Num_Bin = self.spc.shape[1] self.Num_Chn = self.spc.shape[0] - + ''' Se obtiene la constante del RADAR ''' - + self.Pt = Pt self.Gt = Gt self.Gr = Gr @@ -658,30 +676,30 @@ class PrecipitationProc(Operation): self.tauW = tauW self.ThetaT = ThetaT self.ThetaR = ThetaR - + Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR) RadarConstant = 10e-26 * Numerator / Denominator # - + ''' ============================= ''' - - self.spc[0] = (self.spc[0]-dataOut.noise[0]) - self.spc[1] = (self.spc[1]-dataOut.noise[1]) - self.spc[2] = (self.spc[2]-dataOut.noise[2]) - + + self.spc[0] = (self.spc[0]-dataOut.noise[0]) + self.spc[1] = (self.spc[1]-dataOut.noise[1]) + self.spc[2] = (self.spc[2]-dataOut.noise[2]) + self.spc[ numpy.where(self.spc < 0)] = 0 - - SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise)) + + SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise)) SPCmean[ numpy.where(SPCmean < 0)] = 0 - + ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei]) ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei]) ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei]) - + Pr = SPCmean[:,:] - + VelMeteoro = numpy.mean(SPCmean,axis=0) - + D_range = numpy.zeros([self.Num_Bin,self.Num_Hei]) SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei]) N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei]) @@ -690,102 +708,102 @@ class PrecipitationProc(Operation): Z = numpy.zeros(self.Num_Hei) Ze = numpy.zeros(self.Num_Hei) RR = numpy.zeros(self.Num_Hei) - + Range = dataOut.heightList*1000. - + for R in range(self.Num_Hei): - + h = Range[R] + Altitude #Range from ground to radar pulse altitude del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity - + D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3 - + '''NOTA: ETA(n) dn = ETA(f) df - + dn = 1 Diferencial de muestreo df = ETA(n) / ETA(f) - + ''' - + ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA) - + ETAv[:,R]=ETAn[:,R]/dV - + ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R]) - + SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma) - - N_dist[:,R] = ETAn[:,R] / SIGMA[:,R] - + + N_dist[:,R] = ETAn[:,R] / SIGMA[:,R] + DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin]) - + try: popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments) - except: + except: popt01=numpy.zeros(3) popt01[1]= DMoments[1] - + if popt01[1]<0 or popt01[1]>20: popt01[1]=numpy.NaN - - + + V_mean[R]=popt01[1] - + Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18 - + RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate - + Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km) - - - + + + RR2 = (Z/200)**(1/1.6) dBRR = 10*numpy.log10(RR) dBRR2 = 10*numpy.log10(RR2) - + dBZe = 10*numpy.log10(Ze) dBZ = 10*numpy.log10(Z) - + dataOut.data_output = RR[8] dataOut.data_param = numpy.ones([3,self.Num_Hei]) dataOut.channelList = [0,1,2] - + dataOut.data_param[0]=dBZ dataOut.data_param[1]=V_mean dataOut.data_param[2]=RR return dataOut - + def dBZeMODE2(self, dataOut): # Processing for MIRA35C - + NPW = dataOut.NPW COFA = dataOut.COFA - + SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]]) RadarConst = dataOut.RadarConst #frequency = 34.85*10**9 - + ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei])) data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN - + ETA = numpy.sum(SNR,1) - + ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN) - + Ze = numpy.ones([self.Num_Chn, self.Num_Hei] ) - + for r in range(self.Num_Hei): - + Ze[0,r] = ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2) #Ze[1,r] = ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2) - + return Ze - + # def GetRadarConstant(self): -# -# """ +# +# """ # Constants: -# +# # Pt: Transmission Power dB 5kW 5000 # Gt: Transmission Gain dB 24.7 dB 295.1209 # Gr: Reception Gain dB 18.5 dB 70.7945 @@ -794,55 +812,55 @@ class PrecipitationProc(Operation): # tauW: Width of transmission pulse s 4us 4e-6 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087 -# +# # """ -# +# # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR) # RadarConstant = Numerator / Denominator -# +# # return RadarConstant - - - -class FullSpectralAnalysis(Operation): - + + + +class FullSpectralAnalysis(Operation): + """ Function that implements Full Spectral Analysis technique. - - Input: + + Input: self.dataOut.data_pre : SelfSpectra and CrossSpectra data self.dataOut.groupList : Pairlist of channels self.dataOut.ChanDist : Physical distance between receivers - - - Output: - - self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind - - + + + Output: + + self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind + + Parameters affected: Winds, height range, SNR - + """ def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7, minheight=None, maxheight=None): - - self.indice=int(numpy.random.rand()*1000) - + + self.indice=int(numpy.random.rand()*1000) + spc = dataOut.data_pre[0].copy() cspc = dataOut.data_pre[1] - + """Erick: NOTE THE RANGE OF THE PULSE TX MUST BE REMOVED""" SNRspc = spc.copy() SNRspc[:,:,0:7]= numpy.NaN - + """##########################################""" - - + + nChannel = spc.shape[0] nProfiles = spc.shape[1] nHeights = spc.shape[2] - + # first_height = 0.75 #km (ref: data header 20170822) # resolution_height = 0.075 #km ''' @@ -866,37 +884,37 @@ class FullSpectralAnalysis(Operation): ChanDist = dataOut.ChanDist else: ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]]) - + FrecRange = dataOut.spc_range[0] - + data_SNR=numpy.zeros([nProfiles]) noise = dataOut.noise - + dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0] - + dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20 - - + + data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN - + velocityX=[] velocityY=[] - velocityV=[] - + velocityV=[] + dbSNR = 10*numpy.log10(dataOut.data_SNR) dbSNR = numpy.average(dbSNR,0) - + '''***********************************************WIND ESTIMATION**************************************''' - + for Height in range(nHeights): - - if Height >= range_min and Height < range_max: - # error_code unused, yet maybe useful for future analysis. + + if Height >= range_min and Height < range_max: + # error_code unused, yet maybe useful for future analysis. [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit) else: Vzon,Vmer,Vver = 0., 0., numpy.NaN - - + + if abs(Vzon) < 100. and abs(Vzon) > 0. and abs(Vmer) < 100. and abs(Vmer) > 0.: velocityX=numpy.append(velocityX, Vzon) velocityY=numpy.append(velocityY, -Vmer) @@ -904,33 +922,33 @@ class FullSpectralAnalysis(Operation): else: velocityX=numpy.append(velocityX, numpy.NaN) velocityY=numpy.append(velocityY, numpy.NaN) - + if dbSNR[Height] > SNRlimit: velocityV=numpy.append(velocityV, -Vver) # reason for this minus sign -> convention? (taken from Ericks version) else: velocityV=numpy.append(velocityV, numpy.NaN) - - + + '''Change the numpy.array (velocityX) sign when trying to process BLTR data (Erick)''' - data_output[0] = numpy.array(velocityX) - data_output[1] = numpy.array(velocityY) + data_output[0] = numpy.array(velocityX) + data_output[1] = numpy.array(velocityY) data_output[2] = velocityV - - + + dataOut.data_output = data_output - + return dataOut - + def moving_average(self,x, N=2): """ convolution for smoothenig data. note that last N-1 values are convolution with zeroes """ return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):] - + def gaus(self,xSamples,Amp,Mu,Sigma): return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) )) - + def Moments(self, ySamples, xSamples): - '''*** + '''*** Variables corresponding to moments of distribution. Also used as initial coefficients for curve_fit. Vr was corrected. Only a velocity when x is velocity, of course. @@ -939,9 +957,9 @@ class FullSpectralAnalysis(Operation): yNorm = ySamples / Pot x_range = (numpy.max(xSamples)-numpy.min(xSamples)) Vr = numpy.nansum( yNorm * xSamples )*x_range/len(xSamples) # Velocidad radial, mu, corrimiento doppler, primer momento - Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento + Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral - + return numpy.array([Pot, Vr, Desv]) def StopWindEstimation(self, error_code): @@ -954,7 +972,7 @@ class FullSpectralAnalysis(Operation): return Vzon, Vmer, Vver, error_code def AntiAliasing(self, interval, maxstep): - """ + """ function to prevent errors from aliased values when computing phaseslope """ antialiased = numpy.zeros(len(interval))*0.0 @@ -964,8 +982,8 @@ class FullSpectralAnalysis(Operation): for i in range(1,len(antialiased)): - step = interval[i] - interval[i-1] - + step = interval[i] - interval[i-1] + if step > maxstep: copyinterval -= 2*numpy.pi antialiased[i] = copyinterval[i] @@ -973,7 +991,7 @@ class FullSpectralAnalysis(Operation): elif step < maxstep*(-1): copyinterval += 2*numpy.pi antialiased[i] = copyinterval[i] - + else: antialiased[i] = copyinterval[i].copy() @@ -1003,27 +1021,27 @@ class FullSpectralAnalysis(Operation): 3 : SNR to low or velocity to high -> prec. e.g. 4 : at least one Gaussian of cspc exceeds widthlimit 5 : zero out of three cspc Gaussian fits converged - 6 : phase slope fit could not be found + 6 : phase slope fit could not be found 7 : arrays used to fit phase have different length 8 : frequency range is either too short (len <= 5) or very long (> 30% of cspc) """ error_code = 0 - + SPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]]) # for normalized spc values for one height phase = numpy.ones([spc.shape[0],spc.shape[1]]) # phase between channels CSPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) # for normalized cspc values PhaseSlope = numpy.zeros(spc.shape[0]) # slope of the phases, channelwise PhaseInter = numpy.ones(spc.shape[0]) # intercept to the slope of the phases, channelwise - xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range + xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range xVel = AbbsisaRange[2][0:spc.shape[1]] # velocity range SPCav = numpy.average(spc, axis=0)-numpy.average(noise) # spc[0]-noise[0] - + SPCmoments_vel = self.Moments(SPCav, xVel ) # SPCmoments_vel[1] corresponds to vertical velocity and is used to determine if signal corresponds to wind (if .. <3) CSPCmoments = [] - + '''Getting Eij and Nij''' @@ -1038,13 +1056,13 @@ class FullSpectralAnalysis(Operation): spc_norm = spc.copy() # need copy() because untouched spc is needed for normalization of cspc below spc_norm = numpy.where(numpy.isfinite(spc_norm), spc_norm, numpy.NAN) - for i in range(spc.shape[0]): - + for i in range(spc.shape[0]): + spc_sub = spc_norm[i,:] - noise[i] # spc not smoothed here or in previous version. Factor_Norm = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc_sub)) # usually = Freq range / nfft - normalized_spc = spc_sub / (numpy.nansum(numpy.abs(spc_sub)) * Factor_Norm) - + normalized_spc = spc_sub / (numpy.nansum(numpy.abs(spc_sub)) * Factor_Norm) + xSamples = xFrec # the frequency range is taken SPC_Samples[i] = normalized_spc # Normalized SPC values are taken @@ -1055,49 +1073,49 @@ class FullSpectralAnalysis(Operation): only for estimation of width. for normalization of cross spectra, you need initial, unnormalized self-spectra With noise. - Technically, you don't even need to normalize the self-spectra, as you only need the + Technically, you don't even need to normalize the self-spectra, as you only need the width of the peak. However, it was left this way. Note that the normalization has a flaw: due to subtraction of the noise, some values are below zero. Raw "spc" values should be >= 0, as it is the modulus squared of the signals (complex * it's conjugate) """ - SPCMean = numpy.average(SPC_Samples, axis=0) - + SPCMean = numpy.average(SPC_Samples, axis=0) + popt = [1e-10,0,1e-10] SPCMoments = self.Moments(SPCMean, xSamples) - if dbSNR > SNRlimit and numpy.abs(SPCmoments_vel[1]) < 3: + if dbSNR > SNRlimit and numpy.abs(SPCmoments_vel[1]) < 3: try: popt,pcov = curve_fit(self.gaus,xSamples,SPCMean,p0=SPCMoments)#, bounds=(-numpy.inf, [numpy.inf, numpy.inf, 10])). Setting bounds does not make the code faster but only keeps the fit from finding the minimum. if popt[2] > widthlimit: # CONDITION return self.StopWindEstimation(error_code = 1) FitGauss = self.gaus(xSamples,*popt) - + except :#RuntimeError: return self.StopWindEstimation(error_code = 2) else: return self.StopWindEstimation(error_code = 3) - + '''***************************** CSPC Normalization ************************* new section: The Spc spectra are used to normalize the crossspectra. Peaks from precipitation - influence the norm which is not desired. First, a range is identified where the - wind peak is estimated -> sum_wind is sum of those frequencies. Next, the area - around it gets cut off and values replaced by mean determined by the boundary + influence the norm which is not desired. First, a range is identified where the + wind peak is estimated -> sum_wind is sum of those frequencies. Next, the area + around it gets cut off and values replaced by mean determined by the boundary data -> sum_noise (spc is not normalized here, thats why the noise is important) The sums are then added and multiplied by range/datapoints, because you need an integral and not a sum for normalization. - - A norm is found according to Briggs 92. + + A norm is found according to Briggs 92. ''' radarWavelength = 0.6741 # meters - count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa. + count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa. # count_limit_freq = numpy.max(xFrec) channel_integrals = numpy.zeros(3) @@ -1108,11 +1126,11 @@ class FullSpectralAnalysis(Operation): sum over all frequencies in the range around zero Hz @ math.ceil(N_freq/2) ''' N_freq = numpy.count_nonzero(~numpy.isnan(spc[i,:])) - count_limit_int = int(math.ceil( count_limit_freq / numpy.max(xFrec) * (N_freq / 2) )) # gives integer point + count_limit_int = int(math.ceil( count_limit_freq / numpy.max(xFrec) * (N_freq / 2) )) # gives integer point sum_wind = numpy.nansum( spc[i, (math.ceil(N_freq/2) - count_limit_int) : (math.ceil(N_freq / 2) + count_limit_int)] ) #N_freq/2 is where frequency (velocity) is zero, i.e. middle of spectrum. sum_noise = (numpy.mean(spc[i, :4]) + numpy.mean(spc[i, -6:-2]))/2.0 * (N_freq - 2*count_limit_int) channel_integrals[i] = (sum_noise + sum_wind) * (2*numpy.max(xFrec) / N_freq) - + cross_integrals_peak = numpy.zeros(3) # cross_integrals_totalrange = numpy.zeros(3) @@ -1125,45 +1143,45 @@ class FullSpectralAnalysis(Operation): chan_index1 = pairsList[i][1] cross_integrals_peak[i] = channel_integrals[chan_index0]*channel_integrals[chan_index1] - normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_peak[i]) + normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_peak[i]) CSPC_Samples[i] = normalized_cspc ''' Finding cross integrals without subtracting any peaks:''' # FactorNorm0 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index0,:])) # FactorNorm1 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index1,:])) - # cross_integrals_totalrange[i] = (numpy.nansum(spc[chan_index0,:])) * FactorNorm0 * (numpy.nansum(spc[chan_index1,:])) * FactorNorm1 - # normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_totalrange[i]) + # cross_integrals_totalrange[i] = (numpy.nansum(spc[chan_index0,:])) * FactorNorm0 * (numpy.nansum(spc[chan_index1,:])) * FactorNorm1 + # normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_totalrange[i]) # CSPC_Samples[i] = normalized_cspc - - + + phase[i] = numpy.arctan2(CSPC_Samples[i].imag, CSPC_Samples[i].real) CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0]), xSamples), self.Moments(numpy.abs(CSPC_Samples[1]), xSamples), self.Moments(numpy.abs(CSPC_Samples[2]), xSamples)]) - + '''***Sorting out NaN entries***''' CSPCMask01 = numpy.abs(CSPC_Samples[0]) CSPCMask02 = numpy.abs(CSPC_Samples[1]) CSPCMask12 = numpy.abs(CSPC_Samples[2]) - + mask01 = ~numpy.isnan(CSPCMask01) mask02 = ~numpy.isnan(CSPCMask02) mask12 = ~numpy.isnan(CSPCMask12) - + CSPCMask01 = CSPCMask01[mask01] CSPCMask02 = CSPCMask02[mask02] CSPCMask12 = CSPCMask12[mask12] - - popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10] + + popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10] FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0 - + '''*******************************FIT GAUSS CSPC************************************''' - try: + try: popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0]) if popt01[2] > widthlimit: # CONDITION return self.StopWindEstimation(error_code = 4) @@ -1186,53 +1204,53 @@ class FullSpectralAnalysis(Operation): '''************* Getting Fij ***************''' - - #Punto en Eje X de la Gaussiana donde se encuentra el centro -- x-axis point of the gaussian where the center is located - # -> PointGauCenter - GaussCenter = popt[1] + + #Punto en Eje X de la Gaussiana donde se encuentra el centro -- x-axis point of the gaussian where the center is located + # -> PointGauCenter + GaussCenter = popt[1] ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()] PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0] - + #Punto e^-1 hubicado en la Gaussiana -- point where e^-1 is located in the gaussian PeMinus1 = numpy.max(FitGauss) * numpy.exp(-1) FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss" PointFij = numpy.where(FitGauss==FijClosest)[0][0] Fij = numpy.abs(xSamples[PointFij] - xSamples[PointGauCenter]) - + '''********** Taking frequency ranges from mean SPCs **********''' - + #GaussCenter = popt[1] #Primer momento 01 GauWidth = popt[2] * 3/2 #Ancho de banda de Gau01 -- Bandwidth of Gau01 TODO why *3/2? Range = numpy.empty(2) Range[0] = GaussCenter - GauWidth - Range[1] = GaussCenter + GauWidth + Range[1] = GaussCenter + GauWidth #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max) -- Point in x-axis where the bandwidth is located (min:max) ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()] ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()] - + PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0] PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0] - + Range = numpy.array([ PointRangeMin, PointRangeMax ]) - + FrecRange = xFrec[ Range[0] : Range[1] ] - - '''************************** Getting Phase Slope ***************************''' - - for i in range(1,3): # Changed to only compute two - + + '''************************** Getting Phase Slope ***************************''' + + for i in range(1,3): # Changed to only compute two + if len(FrecRange) > 5 and len(FrecRange) < spc.shape[1] * 0.3: # PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=1) #used before to smooth phase with N=3 PhaseRange = phase[i,Range[0]:Range[1]].copy() - + mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange) - + if len(FrecRange) == len(PhaseRange): - - try: + + try: slope, intercept, _, _, _ = stats.linregress(FrecRange[mask], self.AntiAliasing(PhaseRange[mask], 4.5)) PhaseSlope[i] = slope PhaseInter[i] = intercept @@ -1242,49 +1260,49 @@ class FullSpectralAnalysis(Operation): else: return self.StopWindEstimation(error_code = 7) - + else: return self.StopWindEstimation(error_code = 8) - - - + + + '''*** Constants A-H correspond to the convention as in Briggs and Vincent 1992 ***''' '''Getting constant C''' cC=(Fij*numpy.pi)**2 - + '''****** Getting constants F and G ******''' MijEijNij = numpy.array([[Xi02,Eta02], [Xi12,Eta12]]) MijResult0 = (-PhaseSlope[1] * cC) / (2*numpy.pi) - MijResult1 = (-PhaseSlope[2] * cC) / (2*numpy.pi) + MijResult1 = (-PhaseSlope[2] * cC) / (2*numpy.pi) MijResults = numpy.array([MijResult0,MijResult1]) - (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults) - + (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults) + '''****** Getting constants A, B and H ******''' - W01 = numpy.nanmax( FitGauss01 ) - W02 = numpy.nanmax( FitGauss02 ) - W12 = numpy.nanmax( FitGauss12 ) - + W01 = numpy.nanmax( FitGauss01 ) + W02 = numpy.nanmax( FitGauss02 ) + W12 = numpy.nanmax( FitGauss12 ) + WijResult0 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC)) WijResult1 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC)) WijResult2 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC)) - + WijResults = numpy.array([WijResult0, WijResult1, WijResult2]) - - WijEijNij = numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ]) + + WijEijNij = numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ]) (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults) - + VxVy = numpy.array([[cA,cH],[cH,cB]]) VxVyResults = numpy.array([-cF,-cG]) (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults) - + Vzon = Vy Vmer = Vx - + # Vmag=numpy.sqrt(Vzon**2+Vmer**2) # unused # Vang=numpy.arctan2(Vmer,Vzon) # unused - + ''' using frequency as abscissa. Due to three channels, the offzenith angle is zero and Vrad equal to Vver. formula taken from Briggs 92, figure 4. ''' @@ -1295,62 +1313,62 @@ class FullSpectralAnalysis(Operation): error_code = 0 - return Vzon, Vmer, Vver, error_code + return Vzon, Vmer, Vver, error_code class SpectralMoments(Operation): - + ''' Function SpectralMoments() - + Calculates moments (power, mean, standard deviation) and SNR of the signal - + Type of dataIn: Spectra - + Configuration Parameters: - + dirCosx : Cosine director in X axis dirCosy : Cosine director in Y axis - + elevation : azimuth : - + Input: - channelList : simple channel list to select e.g. [2,3,7] + channelList : simple channel list to select e.g. [2,3,7] self.dataOut.data_pre : Spectral data self.dataOut.abscissaList : List of frequencies self.dataOut.noise : Noise level per channel - + Affected: self.dataOut.moments : Parameters per channel self.dataOut.data_SNR : SNR per channel - + ''' - + def run(self, dataOut): - + #dataOut.data_pre = dataOut.data_pre[0] data = dataOut.data_pre[0] absc = dataOut.abscissaList[:-1] noise = dataOut.noise nChannel = data.shape[0] data_param = numpy.zeros((nChannel, 4, data.shape[2])) - + for ind in range(nChannel): data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] ) - + dataOut.moments = data_param[:,1:,:] dataOut.data_SNR = data_param[:,0] dataOut.data_POW = data_param[:,1] dataOut.data_DOP = data_param[:,2] dataOut.data_WIDTH = data_param[:,3] return dataOut - - def __calculateMoments(self, oldspec, oldfreq, n0, + + def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None): - + if (nicoh is None): nicoh = 1 - if (graph is None): graph = 0 + if (graph is None): graph = 0 if (smooth is None): smooth = 0 elif (self.smooth < 3): smooth = 0 @@ -1361,98 +1379,98 @@ class SpectralMoments(Operation): if (aliasing is None): aliasing = 0 if (oldfd is None): oldfd = 0 if (wwauto is None): wwauto = 0 - + if (n0 < 1.e-20): n0 = 1.e-20 - + freq = oldfreq vec_power = numpy.zeros(oldspec.shape[1]) vec_fd = numpy.zeros(oldspec.shape[1]) vec_w = numpy.zeros(oldspec.shape[1]) vec_snr = numpy.zeros(oldspec.shape[1]) - + oldspec = numpy.ma.masked_invalid(oldspec) for ind in range(oldspec.shape[1]): - + spec = oldspec[:,ind] aux = spec*fwindow max_spec = aux.max() m = list(aux).index(max_spec) - - #Smooth + + #Smooth if (smooth == 0): spec2 = spec else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth) - + # Calculo de Momentos bb = spec2[list(range(m,spec2.size))] bb = (bb m): ss1 = m - - valid = numpy.asarray(list(range(int(m + bb0 - ss1 + 1)))) + ss1 + + valid = numpy.asarray(list(range(int(m + bb0 - ss1 + 1)))) + ss1 power = ((spec2[valid] - n0)*fwindow[valid]).sum() fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power) - snr = (spec2.mean()-n0)/n0 - - if (snr < 1.e-20) : + snr = (spec2.mean()-n0)/n0 + + if (snr < 1.e-20) : snr = 1.e-20 - + vec_power[ind] = power vec_fd[ind] = fd vec_w[ind] = w vec_snr[ind] = snr - + moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w)) return moments - + #------------------ Get SA Parameters -------------------------- - + def GetSAParameters(self): #SA en frecuencia 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] - + 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) - + 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 ---------------------------------- - + class SALags(Operation): ''' Function GetMoments() @@ -1465,19 +1483,19 @@ class SALags(Operation): self.dataOut.data_SNR self.dataOut.groupList self.dataOut.nChannels - + Affected: self.dataOut.data_param - + ''' - def run(self, dataOut): + def run(self, dataOut): data_acf = dataOut.data_pre[0] data_ccf = dataOut.data_pre[1] normFactor_acf = dataOut.normFactor[0] normFactor_ccf = dataOut.normFactor[1] pairs_acf = dataOut.groupList[0] pairs_ccf = dataOut.groupList[1] - + nHeights = dataOut.nHeights absc = dataOut.abscissaList noise = dataOut.noise @@ -1488,97 +1506,97 @@ class SALags(Operation): for l in range(len(pairs_acf)): data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:] - + for l in range(len(pairs_ccf)): data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:] - + dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights)) dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc) dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc) return - + # def __getPairsAutoCorr(self, pairsList, nChannels): -# +# # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan -# -# for l in range(len(pairsList)): +# +# for l in range(len(pairsList)): # firstChannel = pairsList[l][0] # secondChannel = pairsList[l][1] -# -# #Obteniendo pares de Autocorrelacion +# +# #Obteniendo pares de Autocorrelacion # if firstChannel == secondChannel: # pairsAutoCorr[firstChannel] = int(l) -# +# # pairsAutoCorr = pairsAutoCorr.astype(int) -# +# # pairsCrossCorr = range(len(pairsList)) # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr) -# +# # return pairsAutoCorr, pairsCrossCorr - + def __calculateTaus(self, data_acf, data_ccf, lagRange): - + lag0 = data_acf.shape[1]/2 #Funcion de Autocorrelacion mean_acf = stats.nanmean(data_acf, axis = 0) - + #Obtencion Indice de TauCross ind_ccf = data_ccf.argmax(axis = 1) #Obtencion Indice de TauAuto ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int') ccf_lag0 = data_ccf[:,lag0,:] - + for i in range(ccf_lag0.shape[0]): ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0) - + #Obtencion de TauCross y TauAuto tau_ccf = lagRange[ind_ccf] tau_acf = lagRange[ind_acf] - + Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0]) - + tau_ccf[Nan1,Nan2] = numpy.nan tau_acf[Nan1,Nan2] = numpy.nan tau = numpy.vstack((tau_ccf,tau_acf)) - + return tau - + def __calculateLag1Phase(self, data, lagTRange): data1 = stats.nanmean(data, axis = 0) lag1 = numpy.where(lagTRange == 0)[0][0] + 1 phase = numpy.angle(data1[lag1,:]) - + return phase - + class SpectralFitting(Operation): ''' Function GetMoments() - + Input: Output: Variables modified: ''' - - def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None): - - + + def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None): + + if path != None: sys.path.append(path) self.dataOut.library = importlib.import_module(file) - + #To be inserted as a parameter groupArray = numpy.array(groupList) -# groupArray = numpy.array([[0,1],[2,3]]) +# groupArray = numpy.array([[0,1],[2,3]]) self.dataOut.groupList = groupArray - + nGroups = groupArray.shape[0] nChannels = self.dataIn.nChannels nHeights=self.dataIn.heightList.size - + #Parameters Array self.dataOut.data_param = None - + #Set constants constants = self.dataOut.library.setConstants(self.dataIn) self.dataOut.constants = constants @@ -1587,24 +1605,24 @@ class SpectralFitting(Operation): ippSeconds = self.dataIn.ippSeconds K = self.dataIn.nIncohInt pairsArray = numpy.array(self.dataIn.pairsList) - + #List of possible combinations listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) indCross = numpy.zeros(len(list(listComb)), dtype = 'int') - + if getSNR: listChannels = groupArray.reshape((groupArray.size)) listChannels.sort() noise = self.dataIn.getNoise() self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels]) - - for i in range(nGroups): + + for i in range(nGroups): coord = groupArray[i,:] - + #Input data array data = self.dataIn.data_spc[coord,:,:]/(M*N) data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) - + #Cross Spectra data array for Covariance Matrixes ind = 0 for pairs in listComb: @@ -1613,9 +1631,9 @@ class SpectralFitting(Operation): ind += 1 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N) dataCross = dataCross**2/K - + for h in range(nHeights): - + #Input d = data[:,h] @@ -1624,7 +1642,7 @@ class SpectralFitting(Operation): ind = 0 for pairs in listComb: #Coordinates in Covariance Matrix - x = pairs[0] + x = pairs[0] y = pairs[1] #Channel Index S12 = dataCross[ind,:,h] @@ -1638,15 +1656,15 @@ class SpectralFitting(Operation): LT=L.T dp = numpy.dot(LT,d) - + #Initial values data_spc = self.dataIn.data_spc[coord,:,h] - + if (h>0)and(error1[3]<5): p0 = self.dataOut.data_param[i,:,h-1] else: p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i)) - + try: #Least Squares minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) @@ -1659,30 +1677,30 @@ class SpectralFitting(Operation): minp = p0*numpy.nan error0 = numpy.nan error1 = p0*numpy.nan - + #Save if self.dataOut.data_param is None: self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan - + self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) self.dataOut.data_param[i,:,h] = minp return - + def __residFunction(self, p, dp, LT, constants): fm = self.dataOut.library.modelFunction(p, constants) fmp=numpy.dot(LT,fm) - + return dp-fmp def __getSNR(self, z, noise): - + avg = numpy.average(z, axis=1) SNR = (avg.T-noise)/noise SNR = SNR.T return SNR - + def __chisq(p,chindex,hindex): #similar to Resid but calculates CHI**2 [LT,d,fm]=setupLTdfm(p,chindex,hindex) @@ -1690,53 +1708,53 @@ class SpectralFitting(Operation): fmp=numpy.dot(LT,fm) chisq=numpy.dot((dp-fmp).T,(dp-fmp)) return chisq - + class WindProfiler(Operation): - + __isConfig = False - + __initime = None __lastdatatime = None __integrationtime = None - + __buffer = None - + __dataReady = False - + __firstdata = None - + n = None - - def __init__(self): + + def __init__(self): Operation.__init__(self) - + def __calculateCosDir(self, elev, azim): zen = (90 - elev)*numpy.pi/180 azim = azim*numpy.pi/180 - cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2))) + cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2))) cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2) - + signX = numpy.sign(numpy.cos(azim)) signY = numpy.sign(numpy.sin(azim)) - + cosDirX = numpy.copysign(cosDirX, signX) cosDirY = numpy.copysign(cosDirY, signY) return cosDirX, cosDirY - + def __calculateAngles(self, theta_x, theta_y, azimuth): - + dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2) zenith_arr = numpy.arccos(dir_cosw) azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180 - + dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr) dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr) - + return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly): - -# + +# if horOnly: A = numpy.c_[dir_cosu,dir_cosv] else: @@ -1750,37 +1768,37 @@ class WindProfiler(Operation): listPhi = phi.tolist() maxid = listPhi.index(max(listPhi)) minid = listPhi.index(min(listPhi)) - - rango = list(range(len(phi))) + + rango = list(range(len(phi))) # rango = numpy.delete(rango,maxid) - + heiRang1 = heiRang*math.cos(phi[maxid]) heiRangAux = heiRang*math.cos(phi[minid]) indOut = (heiRang1 < heiRangAux[0]).nonzero() heiRang1 = numpy.delete(heiRang1,indOut) - + velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) SNR1 = numpy.zeros([len(phi),len(heiRang1)]) - + for i in rango: x = heiRang*math.cos(phi[i]) y1 = velRadial[i,:] f1 = interpolate.interp1d(x,y1,kind = 'cubic') - + x1 = heiRang1 y11 = f1(x1) - + y2 = SNR[i,:] f2 = interpolate.interp1d(x,y2,kind = 'cubic') y21 = f2(x1) - + velRadial1[i,:] = y11 SNR1[i,:] = y21 - + return heiRang1, velRadial1, SNR1 def __calculateVelUVW(self, A, velRadial): - + #Operacion Matricial # velUVW = numpy.zeros((velRadial.shape[1],3)) # for ind in range(velRadial.shape[1]): @@ -1788,27 +1806,27 @@ class WindProfiler(Operation): # velUVW = velUVW.transpose() velUVW = numpy.zeros((A.shape[0],velRadial.shape[1])) velUVW[:,:] = numpy.dot(A,velRadial) - - + + return velUVW - + # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0): - + def techniqueDBS(self, kwargs): """ Function that implements Doppler Beam Swinging (DBS) technique. - + Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, Direction correction (if necessary), Ranges and SNR - + Output: Winds estimation (Zonal, Meridional and Vertical) - + Parameters affected: Winds, height range, SNR """ velRadial0 = kwargs['velRadial'] heiRang = kwargs['heightList'] SNR0 = kwargs['SNR'] - + if 'dirCosx' in kwargs and 'dirCosy' in kwargs: theta_x = numpy.array(kwargs['dirCosx']) theta_y = numpy.array(kwargs['dirCosy']) @@ -1816,7 +1834,7 @@ class WindProfiler(Operation): elev = numpy.array(kwargs['elevation']) azim = numpy.array(kwargs['azimuth']) theta_x, theta_y = self.__calculateCosDir(elev, azim) - azimuth = kwargs['correctAzimuth'] + azimuth = kwargs['correctAzimuth'] if 'horizontalOnly' in kwargs: horizontalOnly = kwargs['horizontalOnly'] else: horizontalOnly = False @@ -1831,22 +1849,22 @@ class WindProfiler(Operation): param = param[arrayChannel,:,:] theta_x = theta_x[arrayChannel] theta_y = theta_y[arrayChannel] - - azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) - heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0) + + azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0) A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly) - + #Calculo de Componentes de la velocidad con DBS winds = self.__calculateVelUVW(A,velRadial1) - + return winds, heiRang1, SNR1 - + def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None): - + nPairs = len(pairs_ccf) posx = numpy.asarray(posx) posy = numpy.asarray(posy) - + #Rotacion Inversa para alinear con el azimuth if azimuth!= None: azimuth = azimuth*math.pi/180 @@ -1855,126 +1873,126 @@ class WindProfiler(Operation): else: posx1 = posx posy1 = posy - + #Calculo de Distancias distx = numpy.zeros(nPairs) disty = numpy.zeros(nPairs) dist = numpy.zeros(nPairs) ang = numpy.zeros(nPairs) - + for i in range(nPairs): distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]] - disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]] + disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]] dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2) ang[i] = numpy.arctan2(disty[i],distx[i]) - + return distx, disty, dist, ang - #Calculo de Matrices + #Calculo de Matrices # nPairs = len(pairs) # ang1 = numpy.zeros((nPairs, 2, 1)) # dist1 = numpy.zeros((nPairs, 2, 1)) -# +# # for j in range(nPairs): # dist1[j,0,0] = dist[pairs[j][0]] # dist1[j,1,0] = dist[pairs[j][1]] # ang1[j,0,0] = ang[pairs[j][0]] # ang1[j,1,0] = ang[pairs[j][1]] -# +# # return distx,disty, dist1,ang1 - + def __calculateVelVer(self, phase, lagTRange, _lambda): Ts = lagTRange[1] - lagTRange[0] velW = -_lambda*phase/(4*math.pi*Ts) - + return velW - + def __calculateVelHorDir(self, dist, tau1, tau2, ang): nPairs = tau1.shape[0] nHeights = tau1.shape[1] - vel = numpy.zeros((nPairs,3,nHeights)) + vel = numpy.zeros((nPairs,3,nHeights)) dist1 = numpy.reshape(dist, (dist.size,1)) - + angCos = numpy.cos(ang) angSin = numpy.sin(ang) - - vel0 = dist1*tau1/(2*tau2**2) + + vel0 = dist1*tau1/(2*tau2**2) vel[:,0,:] = (vel0*angCos).sum(axis = 1) vel[:,1,:] = (vel0*angSin).sum(axis = 1) - + ind = numpy.where(numpy.isinf(vel)) vel[ind] = numpy.nan - + return vel - + # def __getPairsAutoCorr(self, pairsList, nChannels): -# +# # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan -# -# for l in range(len(pairsList)): +# +# for l in range(len(pairsList)): # firstChannel = pairsList[l][0] # secondChannel = pairsList[l][1] -# -# #Obteniendo pares de Autocorrelacion +# +# #Obteniendo pares de Autocorrelacion # if firstChannel == secondChannel: # pairsAutoCorr[firstChannel] = int(l) -# +# # pairsAutoCorr = pairsAutoCorr.astype(int) -# +# # pairsCrossCorr = range(len(pairsList)) # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr) -# +# # return pairsAutoCorr, pairsCrossCorr - + # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor): def techniqueSA(self, kwargs): - - """ + + """ Function that implements Spaced Antenna (SA) technique. - + Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, Direction correction (if necessary), Ranges and SNR - + Output: Winds estimation (Zonal, Meridional and Vertical) - + Parameters affected: Winds """ position_x = kwargs['positionX'] position_y = kwargs['positionY'] azimuth = kwargs['azimuth'] - + if 'correctFactor' in kwargs: correctFactor = kwargs['correctFactor'] else: correctFactor = 1 - + groupList = kwargs['groupList'] pairs_ccf = groupList[1] tau = kwargs['tau'] _lambda = kwargs['_lambda'] - + #Cross Correlation pairs obtained # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels) # pairsArray = numpy.array(pairsList)[pairsCrossCorr] # pairsSelArray = numpy.array(pairsSelected) # pairs = [] -# +# # #Wind estimation pairs obtained # for i in range(pairsSelArray.shape[0]/2): # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0] # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0] # pairs.append((ind1,ind2)) - + indtau = tau.shape[0]/2 tau1 = tau[:indtau,:] tau2 = tau[indtau:-1,:] # tau1 = tau1[pairs,:] # tau2 = tau2[pairs,:] phase1 = tau[-1,:] - + #--------------------------------------------------------------------- - #Metodo Directo + #Metodo Directo distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth) winds = self.__calculateVelHorDir(dist, tau1, tau2, ang) winds = stats.nanmean(winds, axis=0) @@ -1990,97 +2008,97 @@ class WindProfiler(Operation): winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) winds = correctFactor*winds return winds - + def __checkTime(self, currentTime, paramInterval, outputInterval): - + dataTime = currentTime + paramInterval deltaTime = dataTime - self.__initime - + if deltaTime >= outputInterval or deltaTime < 0: self.__dataReady = True - return - + return + def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax): ''' Function that implements winds estimation technique with detected meteors. - + Input: Detected meteors, Minimum meteor quantity to wind estimation - + Output: Winds estimation (Zonal and Meridional) - + Parameters affected: Winds - ''' + ''' #Settings nInt = (heightMax - heightMin)/2 nInt = int(nInt) - winds = numpy.zeros((2,nInt))*numpy.nan - + winds = numpy.zeros((2,nInt))*numpy.nan + #Filter errors error = numpy.where(arrayMeteor[:,-1] == 0)[0] finalMeteor = arrayMeteor[error,:] - + #Meteor Histogram finalHeights = finalMeteor[:,2] hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax)) nMeteorsPerI = hist[0] heightPerI = hist[1] - + #Sort of meteors indSort = finalHeights.argsort() finalMeteor2 = finalMeteor[indSort,:] - + # Calculating winds ind1 = 0 - ind2 = 0 - + ind2 = 0 + for i in range(nInt): nMet = nMeteorsPerI[i] ind1 = ind2 ind2 = ind1 + nMet - + meteorAux = finalMeteor2[ind1:ind2,:] - + if meteorAux.shape[0] >= meteorThresh: vel = meteorAux[:, 6] zen = meteorAux[:, 4]*numpy.pi/180 azim = meteorAux[:, 3]*numpy.pi/180 - + n = numpy.cos(zen) # m = (1 - n**2)/(1 - numpy.tan(azim)**2) # l = m*numpy.tan(azim) l = numpy.sin(zen)*numpy.sin(azim) m = numpy.sin(zen)*numpy.cos(azim) - + A = numpy.vstack((l, m)).transpose() A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose()) windsAux = numpy.dot(A1, vel) - + winds[0,i] = windsAux[0] winds[1,i] = windsAux[1] - + return winds, heightPerI[:-1] - + def techniqueNSM_SA(self, **kwargs): metArray = kwargs['metArray'] heightList = kwargs['heightList'] timeList = kwargs['timeList'] - + rx_location = kwargs['rx_location'] groupList = kwargs['groupList'] azimuth = kwargs['azimuth'] dfactor = kwargs['dfactor'] k = kwargs['k'] - + azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth) d = dist*dfactor #Phase calculation metArray1 = self.__getPhaseSlope(metArray, heightList, timeList) - + metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities - + velEst = numpy.zeros((heightList.size,2))*numpy.nan azimuth1 = azimuth1*numpy.pi/180 - + for i in range(heightList.size): h = heightList[i] indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0] @@ -2093,71 +2111,71 @@ class WindProfiler(Operation): A = numpy.asmatrix(A) A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose() velHor = numpy.dot(A1,velAux) - + velEst[i,:] = numpy.squeeze(velHor) return velEst - + def __getPhaseSlope(self, metArray, heightList, timeList): meteorList = [] #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2 #Putting back together the meteor matrix utctime = metArray[:,0] uniqueTime = numpy.unique(utctime) - + phaseDerThresh = 0.5 ippSeconds = timeList[1] - timeList[0] sec = numpy.where(timeList>1)[0][0] nPairs = metArray.shape[1] - 6 nHeights = len(heightList) - + for t in uniqueTime: metArray1 = metArray[utctime==t,:] # phaseDerThresh = numpy.pi/4 #reducir Phase thresh tmet = metArray1[:,1].astype(int) hmet = metArray1[:,2].astype(int) - + metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1)) metPhase[:,:] = numpy.nan metPhase[:,hmet,tmet] = metArray1[:,6:].T - + #Delete short trails metBool = ~numpy.isnan(metPhase[0,:,:]) heightVect = numpy.sum(metBool, axis = 1) metBool[heightVect phaseDerThresh)) metPhase[phDerAux] = numpy.nan - + #--------------------------METEOR DETECTION ----------------------------------------- indMet = numpy.where(numpy.any(metBool,axis=1))[0] - + for p in numpy.arange(nPairs): phase = metPhase[p,:,:] phDer = metDer[p,:,:] - + for h in indMet: height = heightList[h] phase1 = phase[h,:] #82 phDer1 = phDer[h,:] - + phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap - + indValid = numpy.where(~numpy.isnan(phase1))[0] initMet = indValid[0] endMet = 0 - + for i in range(len(indValid)-1): - + #Time difference inow = indValid[i] inext = indValid[i+1] idiff = inext - inow #Phase difference - phDiff = numpy.abs(phase1[inext] - phase1[inow]) - + phDiff = numpy.abs(phase1[inext] - phase1[inow]) + if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor sizeTrail = inow - initMet + 1 if sizeTrail>3*sec: #Too short meteors @@ -2173,43 +2191,43 @@ class WindProfiler(Operation): vel = slope#*height*1000/(k*d) estAux = numpy.array([utctime,p,height, vel, rsq]) meteorList.append(estAux) - initMet = inext + initMet = inext metArray2 = numpy.array(meteorList) - + return metArray2 - + def __calculateAzimuth1(self, rx_location, pairslist, azimuth0): - + azimuth1 = numpy.zeros(len(pairslist)) dist = numpy.zeros(len(pairslist)) - + for i in range(len(rx_location)): ch0 = pairslist[i][0] ch1 = pairslist[i][1] - + diffX = rx_location[ch0][0] - rx_location[ch1][0] diffY = rx_location[ch0][1] - rx_location[ch1][1] azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi dist[i] = numpy.sqrt(diffX**2 + diffY**2) - + azimuth1 -= azimuth0 return azimuth1, dist - + def techniqueNSM_DBS(self, **kwargs): metArray = kwargs['metArray'] heightList = kwargs['heightList'] - timeList = kwargs['timeList'] + timeList = kwargs['timeList'] azimuth = kwargs['azimuth'] theta_x = numpy.array(kwargs['theta_x']) theta_y = numpy.array(kwargs['theta_y']) - + utctime = metArray[:,0] cmet = metArray[:,1].astype(int) hmet = metArray[:,3].astype(int) SNRmet = metArray[:,4] vmet = metArray[:,5] spcmet = metArray[:,6] - + nChan = numpy.max(cmet) + 1 nHeights = len(heightList) @@ -2225,20 +2243,20 @@ class WindProfiler(Operation): thisH = (h1met>=hmin) & (h1met8) & (vmet<50) & (spcmet<10) indthisH = numpy.where(thisH) - + if numpy.size(indthisH) > 3: - + vel_aux = vmet[thisH] chan_aux = cmet[thisH] cosu_aux = dir_cosu[chan_aux] cosv_aux = dir_cosv[chan_aux] cosw_aux = dir_cosw[chan_aux] - - nch = numpy.size(numpy.unique(chan_aux)) + + nch = numpy.size(numpy.unique(chan_aux)) if nch > 1: A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True) velEst[i,:] = numpy.dot(A,vel_aux) - + return velEst def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs): @@ -2249,39 +2267,39 @@ class WindProfiler(Operation): # noise = dataOut.noise heightList = dataOut.heightList SNR = dataOut.data_SNR - + if technique == 'DBS': - - kwargs['velRadial'] = param[:,1,:] #Radial velocity + + kwargs['velRadial'] = param[:,1,:] #Radial velocity kwargs['heightList'] = heightList kwargs['SNR'] = SNR - + dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.paramInterval - + elif technique == 'SA': - + #Parameters # position_x = kwargs['positionX'] # position_y = kwargs['positionY'] # azimuth = kwargs['azimuth'] -# +# # if kwargs.has_key('crosspairsList'): # pairs = kwargs['crosspairsList'] # else: -# pairs = None -# +# pairs = None +# # if kwargs.has_key('correctFactor'): # correctFactor = kwargs['correctFactor'] # else: # correctFactor = 1 - + # tau = dataOut.data_param # _lambda = dataOut.C/dataOut.frequency # pairsList = dataOut.groupList # nChannels = dataOut.nChannels - + kwargs['groupList'] = dataOut.groupList kwargs['tau'] = dataOut.data_param kwargs['_lambda'] = dataOut.C/dataOut.frequency @@ -2289,30 +2307,30 @@ class WindProfiler(Operation): dataOut.data_output = self.techniqueSA(kwargs) dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.timeInterval - - elif technique == 'Meteors': + + elif technique == 'Meteors': dataOut.flagNoData = True self.__dataReady = False - + if 'nHours' in kwargs: nHours = kwargs['nHours'] - else: + else: nHours = 1 - + if 'meteorsPerBin' in kwargs: meteorThresh = kwargs['meteorsPerBin'] else: meteorThresh = 6 - + if 'hmin' in kwargs: hmin = kwargs['hmin'] else: hmin = 70 if 'hmax' in kwargs: hmax = kwargs['hmax'] else: hmax = 110 - + dataOut.outputInterval = nHours*3600 - + if self.__isConfig == False: # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) #Get Initial LTC time @@ -2320,29 +2338,29 @@ class WindProfiler(Operation): self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() self.__isConfig = True - + if self.__buffer is None: self.__buffer = dataOut.data_param self.__firstdata = copy.copy(dataOut) else: self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) - + self.__checkTime(dataOut.utctime, 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 - + dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) dataOut.flagNoData = False self.__buffer = None - + elif technique == 'Meteors1': dataOut.flagNoData = True self.__dataReady = False - + if 'nMins' in kwargs: nMins = kwargs['nMins'] else: nMins = 20 @@ -2357,7 +2375,7 @@ class WindProfiler(Operation): if 'mode' in kwargs: mode = kwargs['mode'] if 'theta_x' in kwargs: - theta_x = kwargs['theta_x'] + theta_x = kwargs['theta_x'] if 'theta_y' in kwargs: theta_y = kwargs['theta_y'] else: mode = 'SA' @@ -2370,10 +2388,10 @@ class WindProfiler(Operation): freq = 50e6 lamb = C/freq k = 2*numpy.pi/lamb - + timeList = dataOut.abscissaList heightList = dataOut.heightList - + if self.__isConfig == False: dataOut.outputInterval = nMins*60 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) @@ -2384,20 +2402,20 @@ class WindProfiler(Operation): self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() self.__isConfig = True - + if self.__buffer is None: self.__buffer = dataOut.data_param self.__firstdata = copy.copy(dataOut) else: self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) - + self.__checkTime(dataOut.utctime, 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 - + metArray = self.__buffer if mode == 'SA': dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList) @@ -2408,71 +2426,71 @@ class WindProfiler(Operation): self.__buffer = None return - + class EWDriftsEstimation(Operation): - - def __init__(self): - Operation.__init__(self) - + + def __init__(self): + Operation.__init__(self) + def __correctValues(self, heiRang, phi, velRadial, SNR): listPhi = phi.tolist() maxid = listPhi.index(max(listPhi)) minid = listPhi.index(min(listPhi)) - - rango = list(range(len(phi))) + + rango = list(range(len(phi))) # rango = numpy.delete(rango,maxid) - + heiRang1 = heiRang*math.cos(phi[maxid]) heiRangAux = heiRang*math.cos(phi[minid]) indOut = (heiRang1 < heiRangAux[0]).nonzero() heiRang1 = numpy.delete(heiRang1,indOut) - + velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) SNR1 = numpy.zeros([len(phi),len(heiRang1)]) - + for i in rango: x = heiRang*math.cos(phi[i]) y1 = velRadial[i,:] f1 = interpolate.interp1d(x,y1,kind = 'cubic') - + x1 = heiRang1 y11 = f1(x1) - + y2 = SNR[i,:] f2 = interpolate.interp1d(x,y2,kind = 'cubic') y21 = f2(x1) - + velRadial1[i,:] = y11 SNR1[i,:] = y21 - + return heiRang1, velRadial1, SNR1 def run(self, dataOut, zenith, zenithCorrection): heiRang = dataOut.heightList velRadial = dataOut.data_param[:,3,:] SNR = dataOut.data_SNR - + zenith = numpy.array(zenith) - zenith -= zenithCorrection + zenith -= zenithCorrection zenith *= numpy.pi/180 - + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR) - + alp = zenith[0] bet = zenith[1] - + w_w = velRadial1[0,:] w_e = velRadial1[1,:] - - w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp)) - u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp)) - + + w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp)) + u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp)) + winds = numpy.vstack((u,w)) - + dataOut.heightList = heiRang1 dataOut.data_output = winds dataOut.data_SNR = SNR1 - + dataOut.utctimeInit = dataOut.utctime dataOut.outputInterval = dataOut.timeInterval return @@ -2485,11 +2503,11 @@ class NonSpecularMeteorDetection(Operation): data_acf = dataOut.data_pre[0] data_ccf = dataOut.data_pre[1] pairsList = dataOut.groupList[1] - + lamb = dataOut.C/dataOut.frequency tSamp = dataOut.ippSeconds*dataOut.nCohInt paramInterval = dataOut.paramInterval - + nChannels = data_acf.shape[0] nLags = data_acf.shape[1] nProfiles = data_acf.shape[2] @@ -2499,7 +2517,7 @@ class NonSpecularMeteorDetection(Operation): heightList = dataOut.heightList ippSeconds = dataOut.ippSeconds*dataOut.nCohInt*dataOut.nAvg utctime = dataOut.utctime - + dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds) #------------------------ SNR -------------------------------------- @@ -2511,7 +2529,7 @@ class NonSpecularMeteorDetection(Operation): SNR[i] = (power[i]-noise[i])/noise[i] SNRm = numpy.nanmean(SNR, axis = 0) SNRdB = 10*numpy.log10(SNR) - + if mode == 'SA': dataOut.groupList = dataOut.groupList[1] nPairs = data_ccf.shape[0] @@ -2519,22 +2537,22 @@ class NonSpecularMeteorDetection(Operation): phase = numpy.zeros(data_ccf[:,0,:,:].shape) # phase1 = numpy.copy(phase) coh1 = numpy.zeros(data_ccf[:,0,:,:].shape) - + for p in range(nPairs): ch0 = pairsList[p][0] ch1 = pairsList[p][1] ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:]) - phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter -# phase1[p,:,:] = numpy.angle(ccf) #median filter - coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter -# coh1[p,:,:] = numpy.abs(ccf) #median filter + phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter +# phase1[p,:,:] = numpy.angle(ccf) #median filter + coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter +# coh1[p,:,:] = numpy.abs(ccf) #median filter coh = numpy.nanmax(coh1, axis = 0) # struc = numpy.ones((5,1)) # coh = ndimage.morphology.grey_dilation(coh, size=(10,1)) #---------------------- Radial Velocity ---------------------------- phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0) velRad = phaseAux*lamb/(4*numpy.pi*tSamp) - + if allData: boolMetFin = ~numpy.isnan(SNRm) # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0) @@ -2542,31 +2560,31 @@ class NonSpecularMeteorDetection(Operation): #------------------------ Meteor mask --------------------------------- # #SNR mask # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB)) -# +# # #Erase small objects -# boolMet1 = self.__erase_small(boolMet, 2*sec, 5) -# +# boolMet1 = self.__erase_small(boolMet, 2*sec, 5) +# # auxEEJ = numpy.sum(boolMet1,axis=0) # indOver = auxEEJ>nProfiles*0.8 #Use this later # indEEJ = numpy.where(indOver)[0] # indNEEJ = numpy.where(~indOver)[0] -# +# # boolMetFin = boolMet1 -# +# # if indEEJ.size > 0: -# boolMet1[:,indEEJ] = False #Erase heights with EEJ -# +# boolMet1[:,indEEJ] = False #Erase heights with EEJ +# # boolMet2 = coh > cohThresh # boolMet2 = self.__erase_small(boolMet2, 2*sec,5) -# +# # #Final Meteor mask # boolMetFin = boolMet1|boolMet2 - + #Coherence mask boolMet1 = coh > 0.75 struc = numpy.ones((30,1)) boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc) - + #Derivative mask derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0) boolMet2 = derPhase < 0.2 @@ -2583,7 +2601,7 @@ class NonSpecularMeteorDetection(Operation): tmet = coordMet[0] hmet = coordMet[1] - + data_param = numpy.zeros((tmet.size, 6 + nPairs)) data_param[:,0] = utctime data_param[:,1] = tmet @@ -2592,7 +2610,7 @@ class NonSpecularMeteorDetection(Operation): data_param[:,4] = velRad[tmet,hmet] data_param[:,5] = coh[tmet,hmet] data_param[:,6:] = phase[:,tmet,hmet].T - + elif mode == 'DBS': dataOut.groupList = numpy.arange(nChannels) @@ -2600,7 +2618,7 @@ class NonSpecularMeteorDetection(Operation): phase = numpy.angle(data_acf[:,1,:,:]) # phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1)) velRad = phase*lamb/(4*numpy.pi*tSamp) - + #Spectral width # acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1)) # acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1)) @@ -2615,24 +2633,24 @@ class NonSpecularMeteorDetection(Operation): #SNR boolMet1 = (SNRdB>SNRthresh) #SNR mask boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5)) - + #Radial velocity boolMet2 = numpy.abs(velRad) < 20 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5)) - + #Spectral Width boolMet3 = spcWidth < 30 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5)) # boolMetFin = self.__erase_small(boolMet1, 10,5) boolMetFin = boolMet1&boolMet2&boolMet3 - + #Creating data_param coordMet = numpy.where(boolMetFin) cmet = coordMet[0] tmet = coordMet[1] hmet = coordMet[2] - + data_param = numpy.zeros((tmet.size, 7)) data_param[:,0] = utctime data_param[:,1] = cmet @@ -2641,7 +2659,7 @@ class NonSpecularMeteorDetection(Operation): data_param[:,4] = SNR[cmet,tmet,hmet].T data_param[:,5] = velRad[cmet,tmet,hmet].T data_param[:,6] = spcWidth[cmet,tmet,hmet].T - + # self.dataOut.data_param = data_int if len(data_param) == 0: dataOut.flagNoData = True @@ -2651,21 +2669,21 @@ class NonSpecularMeteorDetection(Operation): def __erase_small(self, binArray, threshX, threshY): labarray, numfeat = ndimage.measurements.label(binArray) binArray1 = numpy.copy(binArray) - + for i in range(1,numfeat + 1): auxBin = (labarray==i) auxSize = auxBin.sum() - + x,y = numpy.where(auxBin) widthX = x.max() - x.min() widthY = y.max() - y.min() - + #width X: 3 seg -> 12.5*3 - #width Y: - + #width Y: + if (auxSize < 50) or (widthX < threshX) or (widthY < threshY): binArray1[auxBin] = False - + return binArray1 #--------------- Specular Meteor ---------------- @@ -2675,36 +2693,36 @@ class SMDetection(Operation): Function DetectMeteors() Project developed with paper: HOLDSWORTH ET AL. 2004 - + Input: self.dataOut.data_pre - + centerReceiverIndex: From the channels, which is the center receiver - + hei_ref: Height reference for the Beacon signal extraction tauindex: predefinedPhaseShifts: Predefined phase offset for the voltge signals - + cohDetection: Whether to user Coherent detection or not cohDet_timeStep: Coherent Detection calculation time step cohDet_thresh: Coherent Detection phase threshold to correct phases - + noise_timeStep: Noise calculation time step noise_multiple: Noise multiple to define signal threshold - + multDet_timeLimit: Multiple Detection Removal time limit in seconds multDet_rangeLimit: Multiple Detection Removal range limit in km - + phaseThresh: Maximum phase difference between receiver to be consider a meteor - SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor - + SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor + hmin: Minimum Height of the meteor to use it in the further wind estimations hmax: Maximum Height of the meteor to use it in the further wind estimations azimuth: Azimuth angle correction - + Affected: self.dataOut.data_param - + Rejection Criteria (Errors): 0: No error; analysis OK 1: SNR < SNR threshold @@ -2723,9 +2741,9 @@ class SMDetection(Operation): 14: height ambiguous echo: more then one possible height within 70 to 110 km 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s 16: oscilatory echo, indicating event most likely not an underdense echo - + 17: phase difference in meteor Reestimation - + Data Storage: Meteors for Wind Estimation (8): Utc Time | Range Height @@ -2733,19 +2751,19 @@ class SMDetection(Operation): VelRad errorVelRad Phase0 Phase1 Phase2 Phase3 TypeError - - ''' - + + ''' + def run(self, dataOut, hei_ref = None, tauindex = 0, phaseOffsets = None, - cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, + cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, noise_timeStep = 4, noise_multiple = 4, multDet_timeLimit = 1, multDet_rangeLimit = 3, phaseThresh = 20, SNRThresh = 5, hmin = 50, hmax=150, azimuth = 0, channelPositions = None) : - - + + #Getting Pairslist if channelPositions is None: # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T @@ -2755,53 +2773,53 @@ class SMDetection(Operation): heiRang = dataOut.getHeiRange() #Get Beacon signal - No Beacon signal anymore # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex]) -# +# # if hei_ref != None: # newheis = numpy.where(self.dataOut.heightList>hei_ref) -# - - +# + + #****************REMOVING HARDWARE PHASE DIFFERENCES*************** # see if the user put in pre defined phase shifts voltsPShift = dataOut.data_pre.copy() - + # if predefinedPhaseShifts != None: # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180 -# +# # # 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) -# +# 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 voltsDC = numpy.mean(voltsPShift,1) voltsDC = numpy.mean(voltsDC,1) for i in range(voltsDC.shape[0]): voltsPShift[i] = voltsPShift[i] - voltsDC[i] - - #Don't considerate last heights, theyre used to calculate Hardware Phase Shift + + #Don't considerate last heights, theyre used to calculate Hardware Phase Shift # voltsPShift = voltsPShift[:,:,:newheis[0][0]] - + #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) ********** #Coherent Detection if cohDetection: #use coherent detection to get the net power cohDet_thresh = cohDet_thresh*numpy.pi/180 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh) - + #Non-coherent detection! powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0) #********** END OF COH/NON-COH POWER CALCULATION********************** - + #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS **************** #Get noise noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval) @@ -2811,7 +2829,7 @@ class SMDetection(Operation): #Meteor echoes detection listMeteors = self.__findMeteors(powerNet, signalThresh) #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION ********** - + #************** REMOVE MULTIPLE DETECTIONS (3.5) *************************** #Parameters heiRange = dataOut.getHeiRange() @@ -2821,7 +2839,7 @@ class SMDetection(Operation): #Multiple detection removals listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit) #************ END OF REMOVE MULTIPLE DETECTIONS ********************** - + #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ******************** #Parameters phaseThresh = phaseThresh*numpy.pi/180 @@ -2832,40 +2850,40 @@ class SMDetection(Operation): #Estimation of decay times (Errors N 7, 8, 11) listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency) #******************* END OF METEOR REESTIMATION ******************* - + #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) ************************** #Calculating Radial Velocity (Error N 15) radialStdThresh = 10 - listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval) + listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval) if len(listMeteors4) > 0: #Setting New Array date = dataOut.utctime arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) - + #Correcting phase offset if phaseOffsets != None: phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets) - + #Second Pairslist pairsList = [] pairx = (0,1) pairy = (2,3) pairsList.append(pairx) pairsList.append(pairy) - + jph = numpy.array([0,0,0,0]) h = (hmin,hmax) arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, 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] @@ -2873,73 +2891,73 @@ class SMDetection(Operation): # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax) # error = arrayParameters[:,-1] #********************* END OF PARAMETERS CALCULATION ************************** - - #***************************+ PASS DATA TO NEXT STEP ********************** + + #***************************+ PASS DATA TO NEXT STEP ********************** # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1])) dataOut.data_param = arrayParameters - + if arrayParameters is None: dataOut.flagNoData = True else: dataOut.flagNoData = True - + return - + def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n): - + minIndex = min(newheis[0]) maxIndex = max(newheis[0]) - + voltage = voltage0[:,:,minIndex:maxIndex+1] nLength = voltage.shape[1]/n nMin = 0 nMax = 0 phaseOffset = numpy.zeros((len(pairslist),n)) - + for i in range(n): nMax += nLength phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0])) phaseCCF = numpy.mean(phaseCCF, axis = 2) - phaseOffset[:,i] = phaseCCF.transpose() + phaseOffset[:,i] = phaseCCF.transpose() nMin = nMax # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist) - + #Remove Outliers factor = 2 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5)) dw = numpy.std(wt,axis = 1) dw = dw.reshape((dw.size,1)) - ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor)) + ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor)) phaseOffset[ind] = numpy.nan - phaseOffset = stats.nanmean(phaseOffset, axis=1) - + phaseOffset = stats.nanmean(phaseOffset, axis=1) + return phaseOffset - + def __shiftPhase(self, data, phaseShift): #this will shift the phase of a complex number - dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j) + dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j) return dataShifted - + def __estimatePhaseDifference(self, array, pairslist): nChannel = array.shape[0] nHeights = array.shape[2] numPairs = len(pairslist) # phaseCCF = numpy.zeros((nChannel, 5, nHeights)) phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2])) - + #Correct phases derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:] indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) - - if indDer[0].shape[0] > 0: + + if indDer[0].shape[0] > 0: for i in range(indDer[0].shape[0]): signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]]) phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi - + # for j in range(numSides): # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2]) # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux) -# +# #Linear phaseInt = numpy.zeros((numPairs,1)) angAllCCF = phaseCCF[:,[0,1,3,4],0] @@ -2949,16 +2967,16 @@ class SMDetection(Operation): #Phase Differences phaseDiff = phaseInt - phaseCCF[:,2,:] phaseArrival = phaseInt.reshape(phaseInt.size) - + #Dealias phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival)) # indAlias = numpy.where(phaseArrival > numpy.pi) # phaseArrival[indAlias] -= 2*numpy.pi # indAlias = numpy.where(phaseArrival < -numpy.pi) # phaseArrival[indAlias] += 2*numpy.pi - + return phaseDiff, phaseArrival - + def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh): #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power #find the phase shifts of each channel over 1 second intervals @@ -2968,25 +2986,25 @@ class SMDetection(Operation): numHeights = volts.shape[2] nChannel = volts.shape[0] voltsCohDet = volts.copy() - + pairsarray = numpy.array(pairslist) indSides = pairsarray[:,1] # indSides = numpy.array(range(nChannel)) # indSides = numpy.delete(indSides, indCenter) -# +# # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0) listBlocks = numpy.array_split(volts, numBlocks, 1) - + startInd = 0 endInd = 0 - + for i in range(numBlocks): startInd = endInd - endInd = endInd + listBlocks[i].shape[1] - + endInd = endInd + listBlocks[i].shape[1] + arrayBlock = listBlocks[i] # arrayBlockCenter = listCenter[i] - + #Estimate the Phase Difference phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist) #Phase Difference RMS @@ -2998,21 +3016,21 @@ class SMDetection(Operation): for j in range(indSides.size): arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose()) voltsCohDet[:,startInd:endInd,:] = arrayBlock - + return voltsCohDet - + def __calculateCCF(self, volts, pairslist ,laglist): - + nHeights = volts.shape[2] - nPoints = volts.shape[1] + nPoints = volts.shape[1] voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex') - + for i in range(len(pairslist)): volts1 = volts[pairslist[i][0]] - volts2 = volts[pairslist[i][1]] - + volts2 = volts[pairslist[i][1]] + for t in range(len(laglist)): - idxT = laglist[t] + idxT = laglist[t] if idxT >= 0: vStacked = numpy.vstack((volts2[idxT:,:], numpy.zeros((idxT, nHeights),dtype='complex'))) @@ -3020,10 +3038,10 @@ class SMDetection(Operation): vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'), volts2[:(nPoints + idxT),:])) voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0) - + vStacked = None return voltsCCF - + def __getNoise(self, power, timeSegment, timeInterval): numProfPerBlock = numpy.ceil(timeSegment/timeInterval) numBlocks = int(power.shape[0]/numProfPerBlock) @@ -3032,100 +3050,100 @@ class SMDetection(Operation): listPower = numpy.array_split(power, numBlocks, 0) noise = numpy.zeros((power.shape[0], power.shape[1])) noise1 = numpy.zeros((power.shape[0], power.shape[1])) - + startInd = 0 endInd = 0 - + for i in range(numBlocks): #split por canal startInd = endInd - endInd = endInd + listPower[i].shape[0] - + endInd = endInd + listPower[i].shape[0] + arrayBlock = listPower[i] noiseAux = numpy.mean(arrayBlock, 0) # noiseAux = numpy.median(noiseAux) # noiseAux = numpy.mean(arrayBlock) - noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux - + noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux + noiseAux1 = numpy.mean(arrayBlock) - noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1 - + noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1 + return noise, noise1 - + def __findMeteors(self, power, thresh): nProf = power.shape[0] nHeights = power.shape[1] listMeteors = [] - + for i in range(nHeights): powerAux = power[:,i] threshAux = thresh[:,i] - + indUPthresh = numpy.where(powerAux > threshAux)[0] indDNthresh = numpy.where(powerAux <= threshAux)[0] - + j = 0 - + while (j < indUPthresh.size - 2): if (indUPthresh[j + 2] == indUPthresh[j] + 2): indDNAux = numpy.where(indDNthresh > indUPthresh[j]) indDNthresh = indDNthresh[indDNAux] - + if (indDNthresh.size > 0): indEnd = indDNthresh[0] - 1 indInit = indUPthresh[j] - + meteor = powerAux[indInit:indEnd + 1] indPeak = meteor.argmax() + indInit FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0))) - + listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!! j = numpy.where(indUPthresh == indEnd)[0] + 1 else: j+=1 else: j+=1 - + return listMeteors - + def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit): - - arrayMeteors = numpy.asarray(listMeteors) + + arrayMeteors = numpy.asarray(listMeteors) listMeteors1 = [] - + while arrayMeteors.shape[0] > 0: FLAs = arrayMeteors[:,4] maxFLA = FLAs.argmax() listMeteors1.append(arrayMeteors[maxFLA,:]) - + MeteorInitTime = arrayMeteors[maxFLA,1] MeteorEndTime = arrayMeteors[maxFLA,3] MeteorHeight = arrayMeteors[maxFLA,0] - + #Check neighborhood maxHeightIndex = MeteorHeight + rangeLimit minHeightIndex = MeteorHeight - rangeLimit minTimeIndex = MeteorInitTime - timeLimit maxTimeIndex = MeteorEndTime + timeLimit - + #Check Heights indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex) indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex) indBoth = numpy.where(numpy.logical_and(indTime,indHeight)) - + arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0) - + return listMeteors1 - + def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency): numHeights = volts.shape[2] nChannel = volts.shape[0] - + thresholdPhase = thresh[0] thresholdNoise = thresh[1] thresholdDB = float(thresh[2]) - + thresholdDB1 = 10**(thresholdDB/10) pairsarray = numpy.array(pairslist) indSides = pairsarray[:,1] - + pairslist1 = list(pairslist) pairslist1.append((0,1)) pairslist1.append((3,4)) @@ -3134,31 +3152,31 @@ class SMDetection(Operation): listPowerSeries = [] listVoltageSeries = [] #volts has the war data - + if frequency == 30e6: timeLag = 45*10**-3 else: timeLag = 15*10**-3 lag = numpy.ceil(timeLag/timeInterval) - + for i in range(len(listMeteors)): - + ###################### 3.6 - 3.7 PARAMETERS REESTIMATION ######################### meteorAux = numpy.zeros(16) - + #Loading meteor Data (mHeight, mStart, mPeak, mEnd) mHeight = listMeteors[i][0] mStart = listMeteors[i][1] mPeak = listMeteors[i][2] mEnd = listMeteors[i][3] - + #get the volt data between the start and end times of the meteor meteorVolts = volts[:,mStart:mEnd+1,mHeight] meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1) - + #3.6. Phase Difference estimation phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist) - + #3.7. Phase difference removal & meteor start, peak and end times reestimated #meteorVolts0.- all Channels, all Profiles meteorVolts0 = volts[:,:,mHeight] @@ -3166,15 +3184,15 @@ class SMDetection(Operation): meteorNoise = noise[:,mHeight] meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power - + #Times reestimation mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0] if mStart1.size > 0: mStart1 = mStart1[-1] + 1 - - else: + + else: mStart1 = mPeak - + mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0] if mEndDecayTime1.size == 0: @@ -3182,7 +3200,7 @@ class SMDetection(Operation): else: mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax() - + #meteorVolts1.- all Channels, from start to end meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1] meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1] @@ -3191,17 +3209,17 @@ class SMDetection(Operation): meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1) meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1) ##################### END PARAMETERS REESTIMATION ######################### - + ##################### 3.8 PHASE DIFFERENCE REESTIMATION ######################## # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis - if meteorVolts2.shape[1] > 0: + if meteorVolts2.shape[1] > 0: #Phase Difference re-estimation phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist) meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1]) phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1)) meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting - + #Phase Difference RMS phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1))) powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0) @@ -3216,27 +3234,27 @@ class SMDetection(Operation): #Vectorize meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1] meteorAux[7:11] = phaseDiffint[0:4] - + #Rejection Criterions if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation meteorAux[-1] = 17 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB meteorAux[-1] = 1 - - - else: + + + else: meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd] meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis PowerSeries = 0 - + listMeteors1.append(meteorAux) listPowerSeries.append(PowerSeries) listVoltageSeries.append(meteorVolts1) - - return listMeteors1, listPowerSeries, listVoltageSeries - + + return listMeteors1, listPowerSeries, listVoltageSeries + def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency): - + threshError = 10 #Depending if it is 30 or 50 MHz if frequency == 30e6: @@ -3244,22 +3262,22 @@ class SMDetection(Operation): else: timeLag = 15*10**-3 lag = numpy.ceil(timeLag/timeInterval) - + listMeteors1 = [] - + for i in range(len(listMeteors)): meteorPower = listPower[i] meteorAux = listMeteors[i] - + if meteorAux[-1] == 0: - try: + try: indmax = meteorPower.argmax() indlag = indmax + lag - + y = meteorPower[indlag:] x = numpy.arange(0, y.size)*timeLag - + #first guess a = y[0] tau = timeLag @@ -3268,26 +3286,26 @@ class SMDetection(Operation): y1 = self.__exponential_function(x, *popt) #error estimation error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size)) - + decayTime = popt[1] riseTime = indmax*timeInterval meteorAux[11:13] = [decayTime, error] - + #Table items 7, 8 and 11 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s - meteorAux[-1] = 7 + meteorAux[-1] = 7 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time meteorAux[-1] = 8 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time - meteorAux[-1] = 11 - - + meteorAux[-1] = 11 + + except: - meteorAux[-1] = 11 - - + meteorAux[-1] = 11 + + listMeteors1.append(meteorAux) - + return listMeteors1 #Exponential Function @@ -3295,9 +3313,9 @@ class SMDetection(Operation): def __exponential_function(self, x, a, tau): y = a*numpy.exp(-x/tau) return y - + def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval): - + pairslist1 = list(pairslist) pairslist1.append((0,1)) pairslist1.append((3,4)) @@ -3307,33 +3325,33 @@ class SMDetection(Operation): c = 3e8 lag = numpy.ceil(timeLag/timeInterval) freq = 30e6 - + listMeteors1 = [] - + for i in range(len(listMeteors)): meteorAux = listMeteors[i] if meteorAux[-1] == 0: mStart = listMeteors[i][1] - mPeak = listMeteors[i][2] + mPeak = listMeteors[i][2] mLag = mPeak - mStart + lag - + #get the volt data between the start and end times of the meteor meteorVolts = listVolts[i] meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1) #Get CCF allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2]) - + #Method 2 slopes = numpy.zeros(numPairs) time = numpy.array([-2,-1,1,2])*timeInterval angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0]) - + #Correct phases derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1] indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) - - if indDer[0].shape[0] > 0: + + if indDer[0].shape[0] > 0: for i in range(indDer[0].shape[0]): signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]]) angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi @@ -3342,51 +3360,51 @@ class SMDetection(Operation): for j in range(numPairs): fit = stats.linregress(time, angAllCCF[j,:]) slopes[j] = fit[0] - + #Remove Outlier # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes))) # slopes = numpy.delete(slopes,indOut) # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes))) # slopes = numpy.delete(slopes,indOut) - + radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq) radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq) meteorAux[-2] = radialError meteorAux[-3] = radialVelocity - + #Setting Error #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s - if numpy.abs(radialVelocity) > 200: + if numpy.abs(radialVelocity) > 200: meteorAux[-1] = 15 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity elif radialError > radialStdThresh: meteorAux[-1] = 12 - + listMeteors1.append(meteorAux) return listMeteors1 - + def __setNewArrays(self, listMeteors, date, heiRang): - + #New arrays arrayMeteors = numpy.array(listMeteors) arrayParameters = numpy.zeros((len(listMeteors), 13)) - + #Date inclusion # date = re.findall(r'\((.*?)\)', date) # date = date[0].split(',') # date = map(int, date) -# +# # if len(date)<6: # date.append(0) -# +# # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]] # arrayDate = numpy.tile(date, (len(listMeteors), 1)) arrayDate = numpy.tile(date, (len(listMeteors))) - + #Meteor array # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) - + #Parameters Array arrayParameters[:,0] = arrayDate #Date arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range @@ -3394,13 +3412,13 @@ class SMDetection(Operation): arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases arrayParameters[:,-1] = arrayMeteors[:,-1] #Error - + return arrayParameters - + class CorrectSMPhases(Operation): - + def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None): - + arrayParameters = dataOut.data_param pairsList = [] pairx = (0,1) @@ -3408,49 +3426,49 @@ class CorrectSMPhases(Operation): pairsList.append(pairx) pairsList.append(pairy) jph = numpy.zeros(4) - + phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets) arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets))) - + meteorOps = SMOperations() if channelPositions is None: # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella - + pairslist0, distances = meteorOps.getPhasePairs(channelPositions) h = (hmin,hmax) - + arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph) - + dataOut.data_param = arrayParameters return class SMPhaseCalibration(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, d, phases): gammas = numpy.zeros(2) - + for i in range(len(pairs)): - + pairi = pairs[i] phip3 = phases[:,pairi[0]] @@ -3464,7 +3482,7 @@ class SMPhaseCalibration(Operation): jgamma = numpy.angle(numpy.exp(1j*jgamma)) # jgamma[jgamma>numpy.pi] -= 2*numpy.pi # jgamma[jgamma<-numpy.pi] += 2*numpy.pi - + #Revised distribution jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi)) @@ -3473,39 +3491,39 @@ class SMPhaseCalibration(Operation): 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] - + 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): @@ -3526,16 +3544,16 @@ class SMPhaseCalibration(Operation): 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)): @@ -3543,46 +3561,46 @@ class SMPhaseCalibration(Operation): d2 = d[pairsList[1][1]] d5 = d[pairsList[0][0]] d4 = d[pairsList[0][1]] - + alp2 = alpha_y[iy] #gamma 1 - alp4 = alpha_x[ix] #gamma 0 - + alp4 = alpha_x[ix] #gamma 0 + alp3 = -alp2*d3/d2 - gammas[1] alp5 = -alp4*d5/d4 - gammas[0] # jph[pairy[1]] = alpha_y[iy] -# jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]] - +# jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]] + # jph[pairx[1]] = alpha_x[ix] # jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]] jph[pairsList[0][1]] = alp4 jph[pairsList[0][0]] = alp5 jph[pairsList[1][0]] = alp3 - jph[pairsList[1][1]] = alp2 + jph[pairsList[1][1]] = alp2 jph_array[:,ix,iy] = jph # d = [2.0,2.5,2.5,2.0] - #falta chequear si va a leer bien los meteoros + #falta chequear si va a leer bien los meteoros meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, 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 + phOffset = phOffset*180/numpy.pi return phOffset - - + + def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1): - + dataOut.flagNoData = True - self.__dataReady = False + self.__dataReady = False dataOut.outputInterval = nHours*3600 - + if self.__isConfig == False: # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) #Get Initial LTC time @@ -3590,19 +3608,19 @@ class SMPhaseCalibration(Operation): self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() self.__isConfig = True - + if self.__buffer is None: self.__buffer = dataOut.data_param.copy() else: self.__buffer = numpy.vstack((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 @@ -3624,13 +3642,13 @@ class SMPhaseCalibration(Operation): pairs.append((1,0)) else: pairs.append((0,1)) - + if distances[3] > distances[2]: pairs.append((3,2)) else: pairs.append((2,3)) # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb] - + meteorsArray = self.__buffer error = meteorsArray[:,-1] boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14) @@ -3638,7 +3656,7 @@ class SMPhaseCalibration(Operation): meteorsArray = meteorsArray[ind1,:] meteorsArray[:,-1] = 0 phases = meteorsArray[:,8:12] - + #Calculate Gammas gammas = self.__getGammas(pairs, distances, phases) # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180 @@ -3648,22 +3666,22 @@ class SMPhaseCalibration(Operation): dataOut.data_output = -phasesOff dataOut.flagNoData = False self.__buffer = None - - + + return - + class SMOperations(): - + def __init__(self): - + return - + def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph): - + arrayParameters = arrayParameters0.copy() hmin = h[0] hmax = h[1] - + #Calculate AOA (Error N 3, 4) #JONES ET AL. 1998 AOAthresh = numpy.pi/8 @@ -3671,72 +3689,72 @@ class SMOperations(): phases = -arrayParameters[:,8:12] + jph # phases = numpy.unwrap(phases) arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth) - + #Calculate Heights (Error N 13 and 14) error = arrayParameters[:,-1] Ranges = arrayParameters[:,1] zenith = arrayParameters[:,4] arrayParameters[:,2], 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, directions, error, AOAthresh, azimuth): - + arrayAOA = numpy.zeros((phases.shape[0],3)) cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions) - + 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 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0] error[indError] = 0 #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 + 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 + error[indInvalid] = 4 return arrayAOA, error - + def __getDirectionCosines(self, arrayPhase, pairsList, distances): - + #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): ph0 = arrayPhase[:,pairsList[i][0]] ph1 = arrayPhase[:,pairsList[i][1]] d0 = distances[pairsList[i][0]] d1 = distances[pairsList[i][1]] - - ph0_aux = ph0 + ph1 + + ph0_aux = ph0 + ph1 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux)) # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi -# ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi +# ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi #First Estimation cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1)) - + #Most-Accurate Second Estimation phi1_aux = ph0 - ph1 phi1_aux = phi1_aux.reshape(phi1_aux.size,1) #Direction Cosine 1 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1)) - + #Searching the correct Direction Cosine cosdir0_aux = cosdir0[:,i] cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1) @@ -3745,59 +3763,59 @@ class SMOperations(): 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][:len(ind_h)] 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 + #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0] error[indError] = 0 - indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0] + 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 - + error[indInvalid1] = 13 + return heights, error - + def getPhasePairs(self, channelPositions): chanPos = numpy.array(channelPositions) listOper = list(itertools.combinations(list(range(5)),2)) - + distances = numpy.zeros(4) axisX = [] axisY = [] @@ -3805,15 +3823,15 @@ class SMOperations(): distY = numpy.zeros(3) ix = 0 iy = 0 - + pairX = numpy.zeros((2,2)) pairY = numpy.zeros((2,2)) - + for i in range(len(listOper)): pairi = listOper[i] - + posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:]) - + if posDif[0] == 0: axisY.append(pairi) distY[iy] = posDif[1] @@ -3822,7 +3840,7 @@ class SMOperations(): axisX.append(pairi) distX[ix] = posDif[0] ix += 1 - + for i in range(2): if i==0: dist0 = distX @@ -3830,7 +3848,7 @@ class SMOperations(): else: dist0 = distY axis0 = axisY - + side = numpy.argsort(dist0)[:-1] axis0 = numpy.array(axis0)[side,:] chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0]) @@ -3838,7 +3856,7 @@ class SMOperations(): side = axis1[axis1 != chanC] diff1 = chanPos[chanC,i] - chanPos[side[0],i] diff2 = chanPos[chanC,i] - chanPos[side[1],i] - if diff1<0: + if diff1<0: chan2 = side[0] d2 = numpy.abs(diff1) chan1 = side[1] @@ -3848,7 +3866,7 @@ class SMOperations(): d2 = numpy.abs(diff2) chan1 = side[0] d1 = numpy.abs(diff1) - + if i==0: chanCX = chanC chan1X = chan1 @@ -3860,10 +3878,10 @@ class SMOperations(): chan2Y = chan2 distances[2:4] = numpy.array([d1,d2]) # axisXsides = numpy.reshape(axisX[ix,:],4) -# +# # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0]) # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0]) -# +# # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0] # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0] # channel25X = int(pairX[0,ind25X]) @@ -3872,59 +3890,59 @@ class SMOperations(): # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0] # channel25Y = int(pairY[0,ind25Y]) # channel20Y = int(pairY[1,ind20Y]) - + # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)] - pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)] - + pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)] + return pairslist, distances # 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 +# 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 +# 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 +# phi0_aux[indcsi] -= 2*numpy.pi # indcsi = numpy.where(phi0_aux < -numpy.pi) -# phi0_aux[indcsi] += 2*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) @@ -3933,51 +3951,50 @@ class SMOperations(): # 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] +# +# 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] +# #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 +# error[indInvalid1] = 13 +# +# return heights, error diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index d83d0e0..5bc04b3 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -94,6 +94,12 @@ class SpectraProc(ProcessingUnit): blocksize += dc.size blocksize += spc.size + #print("spc :",spc.shape) + data_wr = None + if self.dataOut.flagWR: + data_wr = fft_volt + blocksize = fft_volt.size + cspc = None pairIndex = 0 if self.dataOut.pairsList != None: @@ -113,16 +119,20 @@ class SpectraProc(ProcessingUnit): pairIndex += 1 blocksize += cspc.size - self.dataOut.data_spc = spc - self.dataOut.data_cspc = cspc - self.dataOut.data_dc = dc - self.dataOut.blockSize = blocksize + self.dataOut.data_spc = spc + self.dataOut.data_cspc = cspc + self.dataOut.data_wr = data_wr + self.dataOut.data_dc = dc + self.dataOut.blockSize = blocksize self.dataOut.flagShiftFFT = False - def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None, shift_fft=False): + def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None, shift_fft=False,flagWR= 0): + + self.dataOut.flagWR = flagWR if self.dataIn.type == "Spectra": self.dataOut.copy(self.dataIn) + if shift_fft: #desplaza a la derecha en el eje 2 determinadas posiciones shift = int(self.dataOut.nFFTPoints/2) @@ -135,7 +145,7 @@ class SpectraProc(ProcessingUnit): return True if self.dataIn.type == "Voltage": - + #print("VOLTAGE INPUT SPECTRA") self.dataOut.flagNoData = True if nFFTPoints == None: @@ -157,6 +167,7 @@ class SpectraProc(ProcessingUnit): nProfiles, self.dataIn.nHeights), dtype='complex') + #print("buffer :",self.buffer.shape) if self.dataIn.flagDataAsBlock: nVoltProfiles = self.dataIn.data.shape[1] @@ -182,6 +193,7 @@ class SpectraProc(ProcessingUnit): self.dataOut.flagNoData = True return 0 else: + #print("Spectra ",self.profIndex) self.buffer[:, self.profIndex, :] = self.dataIn.data.copy() self.profIndex += 1 @@ -191,6 +203,7 @@ class SpectraProc(ProcessingUnit): if self.profIndex == nProfiles: self.__updateSpecFromVoltage() self.__getFft() + #print(" DATAOUT SHAPE SPEC",self.dataOut.data_spc.shape) self.dataOut.flagNoData = False self.firstdatatime = None @@ -291,16 +304,16 @@ class SpectraProc(ProcessingUnit): # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList] self.dataOut.channelList = range(len(channelIndexList)) self.__selectPairsByChannel(channelIndexList) - + return 1 - - + + def selectFFTs(self, minFFT, maxFFT ): """ - Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango + Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango minFFT<= FFT <= maxFFT """ - + if (minFFT > maxFFT): raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT)) @@ -330,20 +343,20 @@ class SpectraProc(ProcessingUnit): self.selectFFTsByIndex(minIndex, maxIndex) return 1 - - + + def setH0(self, h0, deltaHeight = None): - + if not deltaHeight: deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0] - + nHeights = self.dataOut.nHeights - + newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight - + self.dataOut.heightList = newHeiRange - - + + def selectHeights(self, minHei, maxHei): """ Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango @@ -360,7 +373,7 @@ class SpectraProc(ProcessingUnit): 1 si el metodo se ejecuto con exito caso contrario devuelve 0 """ - + if (minHei > maxHei): raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)) @@ -388,7 +401,7 @@ class SpectraProc(ProcessingUnit): maxIndex = len(heights) self.selectHeightsByIndex(minIndex, maxIndex) - + return 1 @@ -436,7 +449,7 @@ class SpectraProc(ProcessingUnit): def selectFFTsByIndex(self, minIndex, maxIndex): """ - + """ if (minIndex < 0) or (minIndex > maxIndex): @@ -459,7 +472,7 @@ class SpectraProc(ProcessingUnit): self.dataOut.data_spc = data_spc self.dataOut.data_cspc = data_cspc self.dataOut.data_dc = data_dc - + self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1]) self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1] self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1] @@ -552,7 +565,7 @@ class SpectraProc(ProcessingUnit): xx_inv = numpy.linalg.inv(xx) xx_aux = xx_inv[0, :] - for ich in range(num_chan): + for ich in range(num_chan): yy = jspectra[ich, ind_vel, :] jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy) @@ -574,12 +587,12 @@ class SpectraProc(ProcessingUnit): return 1 def removeInterference2(self): - + cspc = self.dataOut.data_cspc spc = self.dataOut.data_spc - Heights = numpy.arange(cspc.shape[2]) + Heights = numpy.arange(cspc.shape[2]) realCspc = numpy.abs(cspc) - + for i in range(cspc.shape[0]): LinePower= numpy.sum(realCspc[i], axis=0) Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)] @@ -587,17 +600,17 @@ class SpectraProc(ProcessingUnit): InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 ) InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)] InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)] - - + + InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) ) #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax])) if len(InterferenceRange) 0: + for i in range(self.__nch): + self.__buffer2[i][self.__profIndex2] = data_intensity[i] + self.__bufferV2[i][self.__profIndex2] = data_velocity[i] + #print("Dentro del bucle",self.__buffer2) + self.__profIndex2 += 1 + if self.__profIndex2 == self.m : + data_i = self.__buffer2 + data_v = self.__bufferV2 + #print(data_i.shape) + self.__dataReady = True + self.__profIndex2 = 0 + self.__buffer2 = numpy.zeros((self.__nch,self.m,self.__nHeis)) + self.__bufferV2 = numpy.zeros((self.__nch,self.m,self.__nHeis)) + return data_i,data_v + return None,None + + def pulsePairOp(self,data,datatime=None): + if self.__initime == None: + self.__initime = datatime + + data_intensity,data_velocity = self.pulsePairbyProfiles(data) + self.__lastdatatime = datatime + + if data_intensity is None: + return None,None,None + + avgdatatime = self.__initime + self.__initime = datatime + + return data_intensity,data_velocity,avgdatatime + + def run(self,dataOut,n =None,m=None): + + if not self.isConfig: + self.setup(dataOut = dataOut, n = n, m = m) + self.isConfig = True + + data_intensity,data_velocity,avgdatatime = self.pulsePairOp(dataOut.data_wr,dataOut.utctime) + dataOut.flagNoData = True + + if self.__dataReady: + #print(" DATA " , data_intensity.shape) + #dataOut.data = numpy.array([data_intensity])#aqui amigo revisa + #tmp = numpy.zeros([1,data_intensity.shape[0],data_intensity.shape[1]]) + #tmp[0] = data_intensity + dataOut.data = data_intensity + dataOut.data_velocity = data_velocity + #dataOut.data = tmp + #print(" DATA " , dataOut.data.shape) + dataOut.nIncohInt *= self.n + dataOut.nProfiles = self.m + dataOut.nFFTPoints = self.m + #dataOut.data_intensity = data_intensity + dataOut.PRFbyAngle = self.n + dataOut.utctime = avgdatatime + dataOut.flagNoData = False + #####print("TIEMPO: ",dataOut.utctime) + return dataOut diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 8c75d37..0d769b6 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -8,8 +8,8 @@ from time import time @MPDecorator -class VoltageProc(ProcessingUnit): - +class VoltageProc(ProcessingUnit): + def __init__(self): ProcessingUnit.__init__(self) @@ -115,7 +115,7 @@ class VoltageProc(ProcessingUnit): self.dataOut.data = data # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList] self.dataOut.channelList = range(len(channelIndexList)) - + return 1 def selectHeights(self, minHei=None, maxHei=None): @@ -229,7 +229,7 @@ class VoltageProc(ProcessingUnit): """ Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis] """ - buffer = self.dataOut.data[:, :, 0:int(self.dataOut.nHeights-r)] + buffer = self.dataOut.data[:, :, 0:int(self.dataOut.nHeights-r)] buffer = buffer.reshape(self.dataOut.nChannels, self.dataOut.nProfiles, int(self.dataOut.nHeights/window), window) buffer = numpy.sum(buffer,3) @@ -384,14 +384,16 @@ class CohInt(Operation): """ if not self.__withOverlapping: + print("inside over") self.__buffer += data.copy() self.__profIndex += 1 return #Overlapping data nChannels, nHeis = data.shape + print("show me the light",data.shape) data = numpy.reshape(data, (1, nChannels, nHeis)) - + print(data.shape) #If the buffer is empty then it takes the data value if self.__buffer is None: self.__buffer = data @@ -422,6 +424,7 @@ class CohInt(Operation): """ if not self.__withOverlapping: + #print("ahora que fue") data = self.__buffer n = self.__profIndex @@ -430,6 +433,7 @@ class CohInt(Operation): return data, n + #print("cual funciona") #Integration with Overlapping data = numpy.sum(self.__buffer, axis=0) # print data @@ -445,6 +449,7 @@ class CohInt(Operation): # n = None # print data # raise + #print("beforeputdata") self.putData(data) if self.__profIndex == self.n: @@ -497,8 +502,8 @@ class CohInt(Operation): # print self.__bufferStride[self.__profIndexStride - 1] # raise return self.__bufferStride[self.__profIndexStride - 1] - - + + return None, None def integrate(self, data, datatime=None): @@ -520,7 +525,7 @@ class CohInt(Operation): avgdatatime = self.__initime deltatime = datatime - self.__lastdatatime - + if not self.__withOverlapping: self.__initime = datatime else: @@ -546,7 +551,7 @@ class CohInt(Operation): avgdatatime = (times - 1) * timeInterval + dataOut.utctime self.__dataReady = True return avgdata, avgdatatime - + def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs): if not self.isConfig: @@ -560,12 +565,12 @@ class CohInt(Operation): avgdata, avgdatatime = self.integrateByBlock(dataOut) dataOut.nProfiles /= self.n else: - if stride is None: + if stride is None: avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime) else: avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime) - + # dataOut.timeInterval *= n dataOut.flagNoData = True @@ -670,11 +675,11 @@ class Decoder(Operation): junk = junk.flatten() code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud)) profilesList = range(self.__nProfiles) - - for i in range(self.__nChannels): - for j in profilesList: - self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] - return self.datadecTime + + for i in range(self.__nChannels): + for j in profilesList: + self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] + return self.datadecTime def __convolutionByBlockInFreq(self, data): @@ -691,7 +696,7 @@ class Decoder(Operation): return data - + def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None): if dataOut.flagDecodeData: @@ -722,7 +727,7 @@ class Decoder(Operation): self.__nProfiles = dataOut.nProfiles datadec = None - + if mode == 3: mode = 0 @@ -1105,9 +1110,9 @@ class SplitProfiles(Operation): if shape[2] % n != 0: raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2])) - + new_shape = shape[0], shape[1]*n, int(shape[2]/n) - + dataOut.data = numpy.reshape(dataOut.data, new_shape) dataOut.flagNoData = False @@ -1191,6 +1196,296 @@ class CombineProfiles(Operation): dataOut.ippSeconds *= n return dataOut + + + +class CreateBlockVoltage(Operation): + + isConfig = False + __Index = 0 + bufferShape = None + buffer = None + firstdatatime = None + + def __init__(self,**kwargs): + Operation.__init__(self,**kwargs) + self.isConfig = False + self.__Index = 0 + self.firstdatatime = None + + def setup(self,dataOut, m = None ): + ''' + m= Numero perfiles + ''' + #print("CONFIGURANDO CBV") + self.__nChannels = dataOut.nChannels + self.__nHeis = dataOut.nHeights + shape = dataOut.data.shape #nchannels, nprofiles, nsamples + #print("input nChannels",self.__nChannels) + #print("input nHeis",self.__nHeis) + #print("SETUP CREATE BLOCK VOLTAGE") + #print("input Shape",shape) + #print("dataOut.nProfiles",dataOut.nProfiles) + numberSamples = self.__nHeis + numberProfile = int(m) + dataOut.nProfiles = numberProfile + #print("new numberProfile",numberProfile) + #print("new numberSamples",numberSamples) + + self.bufferShape = shape[0], numberProfile, numberSamples # nchannels,nprofiles,nsamples + self.buffer = numpy.zeros((self.bufferShape)) + self.bufferVel = numpy.zeros((self.bufferShape)) + + def run(self, dataOut, m=None): + #print("RUN") + dataOut.flagNoData = True + dataOut.flagDataAsBlock = False + #print("BLOCK INDEX ",self.__Index) + + if not self.isConfig: + self.setup(dataOut, m= m) + self.isConfig = True + if self.__Index < m: + #print("PROFINDEX BLOCK CBV",self.__Index) + self.buffer[:,self.__Index,:] = dataOut.data + self.bufferVel[:,self.__Index,:] = dataOut.data_velocity + self.__Index += 1 + dataOut.flagNoData = True + + if self.firstdatatime == None: + self.firstdatatime = dataOut.utctime + + if self.__Index == m: + #print("**********************************************") + #print("self.buffer.shape ",self.buffer.shape) + #print("##############",self.firstdatatime) + ##print("*********************************************") + ##print("*********************************************") + ##print("******* nProfiles *******", dataOut.nProfiles) + ##print("*********************************************") + ##print("*********************************************") + dataOut.data = self.buffer + dataOut.data_velocity = self.bufferVel + dataOut.utctime = self.firstdatatime + dataOut.nProfiles = m + self.firstdatatime = None + dataOut.flagNoData = False + dataOut.flagDataAsBlock = True + self.__Index = 0 + dataOut.identifierWR = True + return dataOut + +class PulsePairVoltage(Operation): + ''' + Function PulsePair(Signal Power, Velocity) + The real component of Lag[0] provides Intensity Information + The imag component of Lag[1] Phase provides Velocity Information + + Configuration Parameters: + nPRF = Number of Several PRF + theta = Degree Azimuth angel Boundaries + + Input: + self.dataOut + lag[N] + Affected: + self.dataOut.spc + ''' + isConfig = False + __profIndex = 0 + __initime = None + __lastdatatime = None + __buffer = None + __buffer2 = [] + __buffer3 = None + __dataReady = False + n = None + __nch = 0 + __nHeis = 0 + + def __init__(self,**kwargs): + Operation.__init__(self,**kwargs) + + def setup(self, dataOut, n = None ): + ''' + n= Numero de PRF's de entrada + ''' + self.__initime = None + self.__lastdatatime = 0 + self.__dataReady = False + self.__buffer = 0 + self.__buffer2 = [] + self.__buffer3 = 0 + self.__profIndex = 0 + + self.__nch = dataOut.nChannels + self.__nHeis = dataOut.nHeights + + print("ELVALOR DE n es:", n) + if n == None: + raise ValueError("n should be specified.") + + if n != None: + if n<2: + raise ValueError("n should be greater than 2") + + self.n = n + self.__nProf = n + ''' + if overlapping: + self.__withOverlapping = True + self.__buffer = None + + else: + #print ("estoy sin __withO") + self.__withOverlapping = False + self.__buffer = 0 + self.__buffer2 = [] + self.__buffer3 = 0 + ''' + + def putData(self,data): + ''' + Add a profile to he __buffer and increase in one the __profiel Index + ''' + #print("self.__profIndex :",self.__profIndex) + self.__buffer += data*numpy.conjugate(data) + self.__buffer2.append(numpy.conjugate(data)) + if self.__profIndex > 0: + self.__buffer3 += self.__buffer2[self.__profIndex-1]*data + self.__profIndex += 1 + return + ''' + if not self.__withOverlapping: + #print("Putdata inside over") + self.__buffer += data* numpy.conjugate(data) + self.__buffer2.append(numpy.conjugate(data)) + + if self.__profIndex >0: + self.__buffer3 += self.__buffer2[self.__profIndex-1]*data + self.__profIndex += 1 + return + + if self.__buffer is None: + #print("aqui bro") + self.__buffer = data* numpy.conjugate(data) + self.__buffer2.append(numpy.conjugate(data)) + self.__profIndex += 1 + + return + + if self.__profIndex < self.n: + self.__buffer = numpy.vstack(self.__buffer,data* numpy.conjugate(data)) + self.__buffer2.append(numpy.conjugate(data)) + + if self.__profIndex == 1: + self.__buffer3 = self.__buffer2[self.__profIndex -1] * data + else: + self.__buffer3 = numpy.vstack(self.__buffer3, self.__buffer2[self.profIndex-1]*data) + + self.__profIndex += 1 + return + ''' + + def pushData(self): + ''' + Return the PULSEPAIR and the profiles used in the operation + Affected : self.__profileIndex + ''' + #print("************************************************") + #print("push data int vel n") + data_intensity = self.__buffer/self.n + data_velocity = self.__buffer3/(self.n-1) + n = self.__profIndex + + self.__buffer = 0 + self.__buffer2 = [] + self.__buffer3 = 0 + self.__profIndex = 0 + + return data_intensity, data_velocity,n + ''' + if not self.__withOverlapping: + #print("ahora que fue") + data_intensity = self.__buffer/self.n + data_velocity = self.__buffer3/(self.n-1) + n = self.__profIndex + + self.__buffer = 0 + self.__buffer2 = [] + self.__buffer3 = 0 + self.__profIndex = 0 + return data_intensity, data_velocity,n + + data_intensity = numpy.sum(self.__buffer,axis = 0) + data_velocity = numpy.sum(self.__buffer3,axis = 0) + n = self.__profIndex + #self.__buffer = 0 + #self.__buffer2 = [] + #self.__buffer3 = 0 + #self.__profIndex = 0 + return data_intensity, data_velocity,n + ''' + + def pulsePairbyProfiles(self,data): + + self.__dataReady = False + data_intensity = None + data_velocity = None + #print("beforeputada") + self.putData(data) + #print("ProfileIndex:",self.__profIndex) + if self.__profIndex == self.n: + data_intensity, data_velocity, n = self.pushData() + self.__dataReady = True + #print("-----------------------------------------------") + #print("data_intensity",data_intensity.shape,"data_velocity",data_velocity.shape) + return data_intensity, data_velocity + + def pulsePairOp(self, data, datatime= None): + + if self.__initime == None: + self.__initime = datatime + + data_intensity, data_velocity = self.pulsePairbyProfiles(data) + self.__lastdatatime = datatime + + if data_intensity is None: + return None, None, None + + avgdatatime = self.__initime + deltatime = datatime - self.__lastdatatime + self.__initime = datatime + ''' + if not self.__withOverlapping: + self.__initime = datatime + else: + self.__initime += deltatime + ''' + return data_intensity, data_velocity, avgdatatime + + def run(self, dataOut,n = None, overlapping= False,**kwargs): + + if not self.isConfig: + self.setup(dataOut = dataOut, n = n , **kwargs) + self.isConfig = True + #print("*******************") + #print("print Shape input data:",dataOut.data.shape) + data_intensity, data_velocity, avgdatatime = self.pulsePairOp(dataOut.data, dataOut.utctime) + dataOut.flagNoData = True + + if self.__dataReady: + #print("#------------------------------------------------------") + #print("data_ready",data_intensity.shape) + dataOut.data = data_intensity #valor para plotear RTI + dataOut.nCohInt *= self.n + dataOut.data_intensity = data_intensity #valor para intensidad + dataOut.data_velocity = data_velocity #valor para velocidad + dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo. + dataOut.utctime = avgdatatime + dataOut.flagNoData = False + return dataOut + # import collections # from scipy.stats import mode # diff --git a/schainpy/scripts/USRP_PLOT_THOR.py b/schainpy/scripts/USRP_PLOT_THOR.py index d2e2eb0..4ed49ce 100644 --- a/schainpy/scripts/USRP_PLOT_THOR.py +++ b/schainpy/scripts/USRP_PLOT_THOR.py @@ -25,9 +25,10 @@ controllerObj.setup(id = '191', name='Test_USRP', description=desc) #path = '/media/data/data/vientos/57.2063km/echoes/NCO_Woodman' -path = '/home/soporte/data_hdf5' #### with clock 35.16 db noise - -figpath = '/home/soporte/data_hdf5_imag' +#path = '/home/soporte/data_hdf5' #### with clock 35.16 db noise +path = '/home/alex/WEATHER_DATA/DATA' +figpath = '/home/alex/WEATHER_DATA/DATA/pic' +#figpath = '/home/soporte/data_hdf5_imag' #remotefolder = "/home/wmaster/graficos" ####################################################################### ################# RANGO DE PLOTEO###################################### @@ -95,14 +96,33 @@ procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=rea #opObj11.addParameter(name='nBaud', value='28', format='int') #opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other') -#opObj11.addParameter(name='n', value='100', format='int') +#opObj11.addParameter(name='n', value='10', format='int') + + +opObj11 = procUnitConfObjA.addOperation(name='PulsePair', optype='other') +opObj11.addParameter(name='n', value='10', format='int') + +opObj11 = procUnitConfObjA.addOperation(name='CreateBlockVoltage', optype='other') +opObj11.addParameter(name='m', value='16', format='int') + +procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjA.getId()) +#Not used because the RGB data is obtained directly from the HF Reader. +#opObj21 = procUnitConfObj2.addOperation(name='GetRGBData') + +opObj21 = procUnitConfObj2.addOperation(name='ParamWriter', optype='external') +opObj21.addParameter(name='path', value=figpath+'/NEWData') +opObj21.addParameter(name='blocksPerFile', value='1', format='int') +opObj21.addParameter(name='metadataList',value='heightList',format='list') +opObj21.addParameter(name='dataList',value='data_intensity',format='list') + +''' ####################################################################### ########## OPERACIONES DOMINIO DE LA FRECUENCIA######################## ####################################################################### procUnitConfObjSousySpectra = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId()) -procUnitConfObjSousySpectra.addParameter(name='nFFTPoints', value='100', format='int') -procUnitConfObjSousySpectra.addParameter(name='nProfiles', value='100', format='int') +procUnitConfObjSousySpectra.addParameter(name='nFFTPoints', value='16', format='int') +procUnitConfObjSousySpectra.addParameter(name='nProfiles', value='16', format='int') #procUnitConfObjSousySpectra.addParameter(name='pairsList', value='(0,0),(1,1),(0,1)', format='pairsList') #opObj13 = procUnitConfObjSousySpectra.addOperation(name='removeDC') @@ -174,10 +194,10 @@ opObj11.addParameter(name='save_period', value=10, format='int') #opObj11 = procUnitConfObjSousySpectra.addOperation(name='SpectraWriter', optype='other') #opObj11.addParameter(name='path', value=wr_path) #opObj11.addParameter(name='blocksPerFile', value='50', format='int') +''' print ("Escribiendo el archivo XML") print ("Leyendo el archivo XML") controllerObj.start() -