From c8092aa564c58a1cfc154e85f8e50f151a1f3253 2022-07-07 20:58:38 From: joabAM Date: 2022-07-07 20:58:38 Subject: [PATCH] grafico outliers, grafico de integraciones, funciones de obtención ruído, etc --- diff --git a/schainc/_noise.c b/schainc/_noise.c index 327e99d..a2cc150 100644 --- a/schainc/_noise.c +++ b/schainc/_noise.c @@ -20,7 +20,7 @@ static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) { } double *sortdata = (double*)PyArray_DATA(data_array); int lenOfData = (int)PyArray_SIZE(data_array) ; - double nums_min = lenOfData*0.2; + double nums_min = lenOfData*0.5;//0.2 if (nums_min <= 5) nums_min = 5; double sump = 0; double sumq = 0; diff --git a/schainpy/controller.py b/schainpy/controller.py index 2bf4c02..c73741c 100644 --- a/schainpy/controller.py +++ b/schainpy/controller.py @@ -340,7 +340,7 @@ class Project(Process): idList = list(self.configurations.keys()) id = int(self.id) * 10 - while True: + while 1: id += 1 if str(id) in idList: @@ -633,7 +633,7 @@ class Project(Process): n = len(self.configurations) flag_no_read = False nProc_noRead = 0 - while not err: + while 1: #print("STAR") n_proc = 0 for conf in self.getUnits(): @@ -672,6 +672,7 @@ class Project(Process): if n == 0: err = True + break def run(self): diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 6c1c888..679de9a 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -199,6 +199,7 @@ class JROData(GenericData): codeList = [] azimuthList = [] elevationList = [] + last_noise = None def __str__(self): @@ -641,7 +642,7 @@ class Spectra(JROData): def setValue(self, value): - print("This property should not be initialized") + print("This property should not be initialized", value) return @@ -875,6 +876,7 @@ class Parameters(Spectra): data_param = None # Parameters obtained data_pre = None # Data Pre Parametrization data_SNR = None # Signal to Noise Ratio + data_outlier = None abscissaList = None # Abscissa, can be velocities, lags or time utctimeInit = None # Initial UTC time paramInterval = None # Time interval to calculate Parameters in seconds @@ -889,7 +891,7 @@ class Parameters(Spectra): nAvg = None noise_estimation = None GauSPC = None # Fit gaussian SPC - + max_nIncohInt = 1 def __init__(self): ''' Constructor diff --git a/schainpy/model/graphics/jroplot_base.py b/schainpy/model/graphics/jroplot_base.py index 0f80196..c157aab 100644 --- a/schainpy/model/graphics/jroplot_base.py +++ b/schainpy/model/graphics/jroplot_base.py @@ -580,7 +580,7 @@ class Plot(Operation): self.sender_queue.append(last_time) - while True: + while 1: try: tm = self.sender_queue.popleft() except IndexError: diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index e8bf6da..2510d0f 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -53,7 +53,7 @@ class SnrPlot(RTIPlot): meta = {} data = { - 'snr': 10 * numpy.log10(dataOut.data_snr) + 'snr': 10 * numpy.log10(dataOut.data_snr * dataOut.nIncohInt/ dataOut.max_nIncohInt) } #print(data['snr']) return data, meta @@ -90,7 +90,7 @@ class PowerPlot(RTIPlot): try: data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) except: - pass + pass return data, {} class SpectralWidthPlot(RTIPlot): diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index fec1ba3..78f3838 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -58,7 +58,10 @@ class SpectraPlot(Plot): spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) data['spc'] = spc data['rti'] = dataOut.getPower() - #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter + noise = 10*numpy.log10(dataOut.getNoise()/float(norm)) + data['noise'] = noise[0] + meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) if self.CODE == 'spc_moments': data['moments'] = dataOut.moments @@ -88,8 +91,8 @@ class SpectraPlot(Plot): z = data['spc'] #print(z.shape, x.shape, y.shape) for n, ax in enumerate(self.axes): - #noise = data['noise'][n] - noise=62 + noise = self.data['noise'][n] + #print(noise) if self.CODE == 'spc_moments': mean = data['moments'][n, 1] if ax.firsttime: @@ -268,7 +271,11 @@ class RTIPlot(Plot): data = {} meta = {} data['rti'] = dataOut.getPower() - data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter + noise = 10*numpy.log10(dataOut.getNoise()/float(norm)) + data['noise'] = noise + return data, meta def plot(self): @@ -327,9 +334,7 @@ class RTIPlot(Plot): if self.showprofile: ax.plot_profile.set_data(data[self.CODE][n], self.y) if "noise" in self.data: - - ax.plot_noise.set_data(numpy.repeat( - data['noise'][n], len(self.y)), self.y) + ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y) class CoherencePlot(RTIPlot): @@ -401,12 +406,17 @@ class NoisePlot(Plot): self.titles = ['Noise'] self.colorbar = False self.plots_adjust.update({'right': 0.85 }) + #if not self.titles: + self.titles = ['Noise Plot'] def update(self, dataOut): data = {} meta = {} - noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1) + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter + #noise = 10*numpy.log10(dataOut.getNoise()/norm) + noise = 10*numpy.log10(dataOut.getNoise()) + noise = noise.reshape(dataOut.nChannels, 1) data['noise'] = noise meta['yrange'] = numpy.array([]) @@ -491,11 +501,11 @@ class SpectraCutPlot(Plot): self.nplots = len(self.data.channels) self.ncols = int(numpy.sqrt(self.nplots) + 0.9) self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9) - self.width = 4.2 * self.ncols + 2.5 + self.width = 4.5 * self.ncols + 2.5 self.height = 4.8 * self.nrows self.ylabel = 'Power [dB]' self.colorbar = False - self.plots_adjust.update({'left':0.05, 'hspace':0.3, 'right': 0.9, 'bottom':0.08}) + self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08}) if len(self.selectedHeightsList) > 0: self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight)) @@ -519,7 +529,11 @@ class SpectraCutPlot(Plot): #print(self.height_index) data = {} meta = {} - spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) + + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter + n0 = 10*numpy.log10(dataOut.getNoise()/float(norm)) + + spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - n0 data['spc'] = spc meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) @@ -824,14 +838,15 @@ class NoiselessSpectraPlot(Plot): self.nplots = len(self.data.channels) self.ncols = int(numpy.sqrt(self.nplots) + 0.9) self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9) - self.height = 2.6 * self.nrows + self.height = 3.5 * self.nrows self.cb_label = 'dB' if self.showprofile: - self.width = 4 * self.ncols + self.width = 5.8 * self.ncols else: - self.width = 3.5 * self.ncols - self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08}) + self.width = 4.8* self.ncols + self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12}) + self.ylabel = 'Range [km]' @@ -848,18 +863,13 @@ class NoiselessSpectraPlot(Plot): norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter n0 = 10*numpy.log10(dataOut.getNoise()/float(norm)) - if self.last_noise == None: - self.last_noise = n0 - else: - n0 = (n0*0.2 + self.last_noise*0.8) - self.last_noise = n0 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) data['spc'] = spc - n0 data['rti'] = dataOut.getPower() - n0 - #data['noise'] = noise + # data['noise'] = noise meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) return data, meta @@ -954,18 +964,11 @@ class NoiselessRTIPlot(Plot): norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter n0 = 10*numpy.log10(dataOut.getNoise()/float(norm)) - #print("noise: ",n0, dataOut.normFactor, norm, dataOut.nIncohInt, dataOut.max_nIncohInt) - if self.last_noise == None: - self.last_noise = n0 - else: - n0 = (n0*0.2 + self.last_noise*0.8) - self.last_noise = n0 + data['noise'] = n0[0] data['noiseless_rti'] = dataOut.getPower() - n0 - #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) - #print(noise) return data, meta def plot(self): @@ -976,19 +979,22 @@ class NoiselessRTIPlot(Plot): self.z = numpy.array(self.z, dtype=float) self.z = numpy.ma.masked_invalid(self.z) + try: if self.channelList != None: if len(self.elevationList) > 0 and len(self.azimuthList) > 0: - self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format( + self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format( self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList] else: - self.titles = ['{} Channel {}'.format( + self.titles = ['{} Channel {}'.format( self.CODE.upper(), x) for x in self.channelList] except: if self.channelList.any() != None: self.titles = ['{} Channel {}'.format( self.CODE.upper(), x) for x in self.channelList] + + if self.decimation is None: x, y, z = self.fill_gaps(self.x, self.y, self.z) else: @@ -996,6 +1002,8 @@ class NoiselessRTIPlot(Plot): dummy_var = self.axes #Extrañamente esto actualiza el valor axes #print("plot shapes ", z.shape, x.shape, y.shape) for n, ax in enumerate(self.axes): + + self.zmin = self.zmin if self.zmin else numpy.min(self.z) self.zmax = self.zmax if self.zmax else numpy.max(self.z) data = self.data[-1] @@ -1027,14 +1035,14 @@ class OutliersRTIPlot(Plot): Plot for data_xxxx object ''' - CODE = 'outlier' + CODE = 'outlier_rtc' # Range Time Counts colormap = 'cool' plot_type = 'pcolorbuffer' def setup(self): self.xaxis = 'time' self.ncols = 1 - self.nrows = self.data.shape('outlier')[0] + self.nrows = self.data.shape('outlier_rtc')[0] self.nplots = self.nrows self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94}) @@ -1049,7 +1057,7 @@ class OutliersRTIPlot(Plot): def update(self, dataOut): data = {} - data['outlier'] = dataOut.data_outlier + data['outlier_rtc'] = dataOut.data_outlier meta = {} @@ -1059,7 +1067,7 @@ class OutliersRTIPlot(Plot): # self.data.normalize_heights() self.x = self.data.times self.y = self.data.yrange - self.z = self.data['outlier'] + self.z = self.data['outlier_rtc'] #self.z = numpy.ma.masked_invalid(self.z) @@ -1085,7 +1093,7 @@ class OutliersRTIPlot(Plot): cmap=self.cmaps[n] ) if self.showprofile: - ax.plot_profile = self.pf_axes[n].plot(data['outlier'][n], self.y)[0] + ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0] self.pf_axes[n].set_xlabel('') else: if self.zlimits is not None: @@ -1097,84 +1105,83 @@ class OutliersRTIPlot(Plot): cmap=self.cmaps[n] ) if self.showprofile: - ax.plot_profile.set_data(data['outlier'][n], self.y) + ax.plot_profile.set_data(data['outlier_rtc'][n], self.y) self.pf_axes[n].set_xlabel('') -# class NoiseRTIPlot(Plot): -# ''' -# Plot for data_xxxx object -# ''' -# -# CODE = 'noise' -# colormap = 'jet' -# plot_type = 'pcolorbuffer' -# -# def setup(self): -# self.xaxis = 'time' -# self.ncols = 1 -# self.nrows = self.data.shape('noise')[0] -# self.nplots = self.nrows -# self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94}) -# -# -# if not self.xlabel: -# self.xlabel = 'Time' -# -# self.ylabel = 'Height [km]' -# if not self.titles: -# self.titles = ['noise Ch:{}'.format(x) for x in range(self.nrows)] -# -# def update(self, dataOut): -# -# data = {} -# norm = dataOut.max_nIncohInt*dataOut.nProfiles* dataOut.nCohInt*dataOut.windowOfFilter -# print("max incoh: ",dataOut.max_nIncohInt ) -# n0 = 10*numpy.log10(dataOut.getNoise()/norm) -# data['noise'] = n0 -# -# meta = {} -# -# return data, meta -# -# def plot(self): -# # self.data.normalize_heights() -# self.x = self.data.times -# self.y = self.data.yrange -# self.z = self.data['noise'] -# -# #self.z = numpy.ma.masked_invalid(self.z) -# -# if self.decimation is None: -# x, y, z = self.fill_gaps(self.x, self.y, self.z) -# else: -# x, y, z = self.fill_gaps(*self.decimate()) -# -# for n, ax in enumerate(self.axes): -# -# self.zmax = self.zmax if self.zmax is not None else numpy.max( -# self.z[n]) -# self.zmin = self.zmin if self.zmin is not None else numpy.min( -# self.z[n]) -# data = self.data[-1] -# if ax.firsttime: -# if self.zlimits is not None: -# self.zmin, self.zmax = self.zlimits[n] -# -# ax.plt = ax.pcolormesh(x, y, z[n].T, -# vmin=self.zmin, -# vmax=self.zmax, -# cmap=self.cmaps[n] -# ) -# if self.showprofile: -# ax.plot_profile = self.pf_axes[n].plot(data['noise'][n], self.y)[0] -# else: -# if self.zlimits is not None: -# self.zmin, self.zmax = self.zlimits[n] -# ax.collections.remove(ax.collections[0]) -# ax.plt = ax.pcolormesh(x, y, z[n].T , -# vmin=self.zmin, -# vmax=self.zmax, -# cmap=self.cmaps[n] -# ) -# if self.showprofile: -# ax.plot_profile.set_data(data['noise'][n], self.y) +class NIncohIntRTIPlot(Plot): + ''' + Plot for data_xxxx object + ''' + + CODE = 'integrations_rtc' # Range Time Counts + colormap = 'BuGn' + plot_type = 'pcolorbuffer' + + def setup(self): + self.xaxis = 'time' + self.ncols = 1 + self.nrows = self.data.shape('integrations_rtc')[0] + self.nplots = self.nrows + self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94}) + + + if not self.xlabel: + self.xlabel = 'Time' + + self.ylabel = 'Height [km]' + if not self.titles: + self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)] + + def update(self, dataOut): + + data = {} + data['integrations_rtc'] = dataOut.nIncohInt + + meta = {} + + return data, meta + + def plot(self): + # self.data.normalize_heights() + self.x = self.data.times + self.y = self.data.yrange + self.z = self.data['integrations_rtc'] + + #self.z = numpy.ma.masked_invalid(self.z) + + if self.decimation is None: + x, y, z = self.fill_gaps(self.x, self.y, self.z) + else: + x, y, z = self.fill_gaps(*self.decimate()) + + for n, ax in enumerate(self.axes): + + self.zmax = self.zmax if self.zmax is not None else numpy.max( + self.z[n]) + self.zmin = self.zmin if self.zmin is not None else numpy.min( + self.z[n]) + data = self.data[-1] + if ax.firsttime: + if self.zlimits is not None: + self.zmin, self.zmax = self.zlimits[n] + + ax.plt = ax.pcolormesh(x, y, z[n].T, + vmin=self.zmin, + vmax=self.zmax, + cmap=self.cmaps[n] + ) + if self.showprofile: + ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0] + self.pf_axes[n].set_xlabel('') + else: + if self.zlimits is not None: + self.zmin, self.zmax = self.zlimits[n] + ax.collections.remove(ax.collections[0]) + ax.plt = ax.pcolormesh(x, y, z[n].T , + vmin=self.zmin, + vmax=self.zmax, + cmap=self.cmaps[n] + ) + if self.showprofile: + ax.plot_profile.set_data(data['integrations_rtc'][n], self.y) + self.pf_axes[n].set_xlabel('') diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index d3db918..6679df9 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -465,7 +465,8 @@ class HDFWriter(Operation): def run(self, dataOut,**kwargs): - self.dataOut = dataOut.copy() + self.dataOut = dataOut + if not(self.isConfig): self.setup(**kwargs) diff --git a/schainpy/model/proc/jroproc_base.py b/schainpy/model/proc/jroproc_base.py index 4390e56..87b8b67 100644 --- a/schainpy/model/proc/jroproc_base.py +++ b/schainpy/model/proc/jroproc_base.py @@ -59,7 +59,7 @@ class ProcessingUnit(object): if mybool: - #print("run jeje") + #print("run yeah") self.run(**kwargs) else: if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error: @@ -210,7 +210,7 @@ def MPDecorator(BaseClass): def run(self): - while True: + while 1: dataOut = self.queue.get() diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 4d8f0af..b63112d 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -101,7 +101,7 @@ class ParametersProc(ProcessingUnit): def run(self): # print("run parameter proc") #---------------------- Voltage Data --------------------------- - + #print(self.dataIn.flagNoData) if self.dataIn.type == "Voltage": self.__updateObjFromInput() @@ -132,10 +132,12 @@ class ParametersProc(ProcessingUnit): self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc] self.dataOut.data_spc = self.dataIn.data_spc self.dataOut.data_cspc = self.dataIn.data_cspc + self.dataOut.data_outlier = self.dataIn.data_outlier self.dataOut.nProfiles = self.dataIn.nProfiles self.dataOut.nIncohInt = self.dataIn.nIncohInt self.dataOut.nFFTPoints = self.dataIn.nFFTPoints self.dataOut.ippFactor = self.dataIn.ippFactor + self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() self.dataOut.ipp = self.dataIn.ipp self.dataOut.abscissaList = self.dataIn.getVelRange(1) @@ -1472,9 +1474,9 @@ class SpectralFitting(Operation): def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None): - print("run SpectralFitting") - self.dataOut = dataOut.copy() - print(self.dataOut.nIncohInt) + #print("run SpectralFitting") + self.dataOut = dataOut#.copy() + #print(self.dataOut.nIncohInt) if path != None: sys.path.append(path) self.dataOut.library = importlib.import_module(file) diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index e4a205b..b9e31e2 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -21,7 +21,7 @@ from schainpy.model.data import _noise from schainpy.utils import log import matplotlib.pyplot as plt -from scipy.optimize import curve_fit +#from scipy.optimize import curve_fit class SpectraProc(ProcessingUnit): @@ -126,6 +126,11 @@ class SpectraProc(ProcessingUnit): def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False): #print("run spc proc") + try: + type = self.dataIn.type.decode("utf-8") + self.dataIn.type = type + except: + pass if self.dataIn.type == "Spectra": try: @@ -209,9 +214,42 @@ class SpectraProc(ProcessingUnit): self.firstdatatime = None self.profIndex = 0 + elif self.dataIn.type == "Parameters": + + self.dataOut.data_spc = self.dataIn.data_spc + self.dataOut.data_cspc = self.dataIn.data_cspc + self.dataOut.data_outlier = self.dataIn.data_outlier + self.dataOut.nProfiles = self.dataIn.nProfiles + self.dataOut.nIncohInt = self.dataIn.nIncohInt + self.dataOut.nFFTPoints = self.dataIn.nFFTPoints + self.dataOut.ippFactor = self.dataIn.ippFactor + self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt + self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() + self.dataOut.ipp = self.dataIn.ipp + #self.dataOut.abscissaList = self.dataIn.getVelRange(1) + #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 + if hasattr(self.dataIn, 'channelList'): + self.dataOut.channelList = self.dataIn.channelList + if hasattr(self.dataIn, 'pairsList'): + self.dataOut.pairsList = self.dataIn.pairsList + self.dataOut.groupList = self.dataIn.pairsList + + self.dataOut.flagNoData = False + + if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels + self.dataOut.ChanDist = self.dataIn.ChanDist + else: self.dataOut.ChanDist = None + + #if hasattr(self.dataIn, 'VelRange'): #Velocities range + # self.dataOut.VelRange = self.dataIn.VelRange + #else: self.dataOut.VelRange = None + + else: - raise ValueError("The type of input object '%s' is not valid".format( + raise ValueError("The type of input object {} is not valid".format( self.dataIn.type)) @@ -497,14 +535,16 @@ class removeDC(Operation): return self.dataOut -class getNoise(Operation): - warnings = False +class getNoiseB(Operation): + + __slots__ =('offset','warnings', 'isConfig', 'minIndex','maxIndex','minIndexFFT','maxIndexFFT') def __init__(self): Operation.__init__(self) + self.isConfig = False + + def setup(self, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False): - def run(self, dataOut, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False): - self.dataOut = dataOut self.warnings = warnings if minHei == None: minHei = self.dataOut.heightList[0] @@ -620,24 +660,111 @@ class getNoise(Operation): except: maxIndexFFT = len( self.dataOut.getFreqRange(1)) + self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT + self.isConfig = True + if offset!=None: + self.offset = 10**(offset/10) + #print("config getNoise Done") + + def run(self, dataOut, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False): + self.dataOut = dataOut + + if not self.isConfig: + self.setup(offset, minHei, maxHei,minVel, maxVel, minFreq, maxFreq, warnings) self.dataOut.noise_estimation = None noise = None if self.dataOut.type == 'Voltage': - noise = self.dataOut.getNoise(ymin_index=minIndex, ymax_index=maxIndex) + noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex) #print(minIndex, maxIndex,minIndexVel, maxIndexVel) elif self.dataOut.type == 'Spectra': - noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex) + + noise = numpy.zeros( self.dataOut.nChannels) + for channel in range( self.dataOut.nChannels): + norm = self.dataOut.max_nIncohInt/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex] + #print("norm nIncoh: ", norm ) + daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex] + daux = numpy.multiply(daux, norm) + #print("offset: ", self.offset, 10*numpy.log10(self.offset)) + #noise[channel] = self.getNoiseByMean(daux)/self.offset + noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset + + #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex) else: - noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex) + noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex) self.dataOut.noise_estimation = noise.copy() # dataOut.noise #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64)) + #print(self.dataOut.flagNoData) return self.dataOut + def getNoiseByMean(self,data): + #data debe estar ordenado + data = numpy.mean(data,axis=1) + sortdata = numpy.sort(data, axis=None) + #sortID=data.argsort() + #print(data.shape) + + pnoise = None + j = 0 + + mean = numpy.mean(sortdata) + min = numpy.min(sortdata) + delta = mean - min + indexes = numpy.where(sortdata > (mean+delta))[0] #only array of indexes + #print(len(indexes)) + if len(indexes)==0: + pnoise = numpy.mean(sortdata) + else: + j = indexes[0] + pnoise = numpy.mean(sortdata[0:j]) + + # from matplotlib import pyplot as plt + # plt.plot(sortdata) + # plt.vlines(j,(pnoise-delta),(pnoise+delta), color='r') + # plt.show() + #print("noise: ", 10*numpy.log10(pnoise)) + return pnoise + + def getNoiseByHS(self,data, navg): + #data debe estar ordenado + #data = numpy.mean(data,axis=1) + sortdata = numpy.sort(data, axis=None) + + lenOfData = len(sortdata) + nums_min = lenOfData*0.05 + + if nums_min <= 5: + + nums_min = 5 + + sump = 0. + sumq = 0. + + j = 0 + cont = 1 + + while((cont == 1)and(j < lenOfData)): + + sump += sortdata[j] + sumq += sortdata[j]**2 + #sumq -= sump**2 + if j > nums_min: + rtest = float(j)/(j-1) + 1.0/0.1 + #if ((sumq*j) > (sump**2)): + if ((sumq*j) > (rtest*sump**2)): + j = j - 1 + sump = sump - sortdata[j] + sumq = sumq - sortdata[j]**2 + cont = 0 + + j += 1 + + lnoise = sump / j + + return lnoise -# import matplotlib.pyplot as plt def fit_func( x, a0, a1, a2): #, a3, a4, a5): z = (x - a1) / a2 @@ -779,7 +906,7 @@ class CleanRayleigh(Operation): #index = tini.tm_hour*12+tini.tm_min/5 ''' - REVISAR + #REVISAR ''' # jspc = jspc/self.nFFTPoints/self.normFactor # jcspc = jcspc/self.nFFTPoints/self.normFactor @@ -1396,7 +1523,7 @@ class IntegrationFaradaySpectra(Operation): data_dc = self.__buffer_dc #(CH, HEIGH) self.maxProfilesInt = self.__profIndex - n = self.__profIndex - self.dataOutliers + n = self.__profIndex - self.dataOutliers # n becomes a matrix self.__buffer_spc = [] self.__buffer_cspc = [] @@ -1455,11 +1582,11 @@ class IntegrationFaradaySpectra(Operation): return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75): - self.dataOut = dataOut.copy() + self.dataOut = dataOut if n == 1: return self.dataOut - self.dataOut.flagNoData = True + if self.dataOut.nChannels == 1: self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS #print(self.dataOut.data_spc.shape, self.dataOut.data_cspc) @@ -1483,7 +1610,7 @@ class IntegrationFaradaySpectra(Operation): self.dataOut.dataLag_spc, self.dataOut.dataLag_cspc, self.dataOut.dataLag_dc) - + self.dataOut.flagNoData = True if self.__dataReady: if not self.ByLags: @@ -1513,6 +1640,7 @@ class IntegrationFaradaySpectra(Operation): self.dataOut.utctime = avgdatatime self.dataOut.flagNoData = False #print("Faraday Integration DONE...") + #print(self.dataOut.flagNoData) return self.dataOut class removeInterference(Operation): @@ -1896,7 +2024,7 @@ class IncohInt(Operation): if dataOut.flagNoData == True: return dataOut - + dataOut.flagNoData = True if not self.isConfig: @@ -1910,7 +2038,7 @@ class IncohInt(Operation): self.incohInt += dataOut.nIncohInt self.nOutliers += dataOut.data_outlier if self.__dataReady: - + #print("prof: ",dataOut.max_nIncohInt,self.__profIndex) dataOut.data_spc = avgdata_spc dataOut.data_cspc = avgdata_cspc dataOut.data_dc = avgdata_dc @@ -1918,7 +2046,7 @@ class IncohInt(Operation): dataOut.data_outlier = self.nOutliers dataOut.utctime = avgdatatime dataOut.flagNoData = False - dataOut.max_nIncohInt = self.__profIndex + dataOut.max_nIncohInt += self.__profIndex self.incohInt = 0 self.nOutliers = 0 self.__profIndex = 0 diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 2d59d67..3acc37c 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -6,7 +6,7 @@ from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon from schainpy.utils import log from schainpy.model.io.utils import getHei_index from time import time -import datetime +#import datetime import numpy #import copy from schainpy.model.data import _noise @@ -1742,6 +1742,7 @@ class SSheightProfiles(Operation): dataOut.step = self.step #print(numpy.shape(dataOut.data)) #exit(1) + #print("new data shape and time:", dataOut.data.shape, dataOut.utctime) return dataOut ################################################################################3############################3 @@ -1754,14 +1755,17 @@ class SSheightProfiles2(Operation): Procesa por perfiles y por bloques ''' - step = None - nsamples = None + bufferShape = None profileShape = None sshProfiles = None profileIndex = None - deltaHeight = None - init_range = None + #nsamples = None + #step = None + #deltaHeight = None + #init_range = None + __slots__ = ('step', 'nsamples', 'deltaHeight', 'init_range', 'isConfig', '__nChannels', + '__nProfiles', '__nHeis', 'deltaHeight', 'new_nHeights') def __init__(self, **kwargs): @@ -1786,12 +1790,12 @@ class SSheightProfiles2(Operation): self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] self.init_range = dataOut.heightList[0] #numberProfile = self.nsamples - self.numberSamples = (self.__nHeis - self.nsamples)/self.step + numberSamples = (self.__nHeis - self.nsamples)/self.step - self.new_nHeights = self.numberSamples + self.new_nHeights = numberSamples - self.bufferShape = int(self.__nChannels), int(self.numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles - self.profileShape = int(self.__nChannels), int(self.nsamples), int(self.numberSamples) # nchannels, nprofiles, nsamples + self.bufferShape = int(self.__nChannels), int(numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles + self.profileShape = int(self.__nChannels), int(self.nsamples), int(numberSamples) # nchannels, nprofiles, nsamples self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex) self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex) @@ -1878,7 +1882,8 @@ class SSheightProfiles2(Operation): dataOut.flagDataAsBlock = True dataBlock = None - #print("new data shape:", dataOut.data.shape) + + #print("new data shape:", dataOut.data.shape, dataOut.utctime) return dataOut @@ -1893,27 +1898,28 @@ class removeProfileByFaradayHS(Operation): ''' ''' - isConfig = False - n = None - - __profIndex = 0 + #isConfig = False + #n = None - __dataReady = False + #__dataReady = False __buffer_data = [] __buffer_times = [] - __initime = None - __count_exec = 0 - __profIndex = 0 + #__initime = None + #__count_exec = 0 + #__profIndex = 0 buffer = None - lenProfileOut = 1 + #lenProfileOut = 1 + + #init_prof = 0 + #end_prof = 0 - init_prof = 0 - end_prof = 0 - n_profiles = 0 - first_utcBlock = None + #first_utcBlock = None outliers_IDs_list = [] - __dh = 0 + #__dh = 0 + __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights', + '__dh','first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels', + '__count_exec','__initime','__dataReady','__ipp') def __init__(self, **kwargs): Operation.__init__(self, **kwargs) @@ -1932,12 +1938,13 @@ class removeProfileByFaradayHS(Operation): self.thHistOutlier = thHistOutlier self.__profIndex = 0 self.buffer = None - self.lenProfileOut = 1 + self._ipp = dataOut.ippSeconds self.n_prof_released = 0 self.heightList = dataOut.heightList self.init_prof = 0 self.end_prof = 0 - self.n_profiles = 0 + self.__count_exec = 0 + self.__profIndex = 0 self.first_utcBlock = None self.__dh = dataOut.heightList[1] - dataOut.heightList[0] minHei = minHei @@ -1948,6 +1955,8 @@ class removeProfileByFaradayHS(Operation): maxHei = dataOut.heightList[-1] self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList) + self.nChannels = dataOut.nChannels + self.nHeights = dataOut.nHeights def filterSatsProfiles(self): data = self.__buffer_data @@ -2124,7 +2133,7 @@ class removeProfileByFaradayHS(Operation): def releaseBlock(self): if self.n % self.lenProfileOut != 0: - raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles)) + raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n)) return None data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt @@ -2139,9 +2148,9 @@ class removeProfileByFaradayHS(Operation): return data def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,th_hist_outlier=3,minHei=None, maxHei=None): - #print("run op buffer 2D") - self.nChannels = dataOut.nChannels - self.nHeights = dataOut.nHeights + #print("run op buffer 2D",dataOut.ippSeconds) + # self.nChannels = dataOut.nChannels + # self.nHeights = dataOut.nHeights if not self.isConfig: #print("init p idx: ", dataOut.profileIndex ) @@ -2158,9 +2167,7 @@ class removeProfileByFaradayHS(Operation): self.lenProfileOut = nProfilesOut dataOut.flagNoData = False #print("tp 2 ",dataOut.data.shape) - #(ch, self.n_profiles, nh) = self.buffer.shape - #print("tp 3 ",self.dataOut.data.shape) - #print("rel: ",self.buffer[:,-1,:]) + self.init_prof = 0 self.end_prof = self.lenProfileOut @@ -2174,12 +2181,13 @@ class removeProfileByFaradayHS(Operation): if numpy.isin(self.n_prof_released, self.outliers_IDs_list): #print("omitting: ", self.n_prof_released) dataOut.flagNoData = True - - dataOut.utctime = self.first_utcBlock + self.init_prof*dataOut.ippSeconds + dataOut.ippSeconds = self._ipp + dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp + # print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds) #dataOut.data = self.releaseBlock() #########################################################3 if self.n % self.lenProfileOut != 0: - raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles)) + raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n)) return None dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt