diff --git a/schainc/_noise.c b/schainc/_noise.c index 38fedc9..327e99d 100644 --- a/schainc/_noise.c +++ b/schainc/_noise.c @@ -66,7 +66,7 @@ static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) { } double *sortdata = (double*)PyArray_DATA(data_array); int lenOfData = (int)PyArray_SIZE(data_array) ; - double nums_min = lenOfData*0.75; + double nums_min = lenOfData*0.85; if (nums_min <= 5) nums_min = 5; double sump = 0; double sumq = 0; diff --git a/schainpy/controller.py b/schainpy/controller.py index 9cbf8fe..2bf4c02 100644 --- a/schainpy/controller.py +++ b/schainpy/controller.py @@ -661,12 +661,10 @@ class Project(Process): elif ok == 'no_Read' and (not flag_no_read): nProc_noRead = n_proc - 1 flag_no_read = True - print("No read") continue elif ok == 'new_Read': nProc_noRead = 0 flag_no_read = False - print("read again") continue elif not ok: #print("not ok",ok) diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index fe08eae..6c1c888 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -78,7 +78,7 @@ def hildebrand_sekhon(data, navg): sortdata = numpy.sort(data, axis=None) ''' lenOfData = len(sortdata) - nums_min = lenOfData*0.2 + nums_min = lenOfData*0.5 if nums_min <= 5: @@ -106,6 +106,7 @@ def hildebrand_sekhon(data, navg): j += 1 lnoise = sump / j + return lnoise ''' return _noise.hildebrand_sekhon(sortdata, navg) @@ -380,7 +381,7 @@ class Voltage(JROData): self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt', 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp'] - def getNoisebyHildebrand(self, channel=None): + def getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None): """ Determino el nivel de ruido usando el metodo Hildebrand-Sekhon @@ -389,10 +390,10 @@ class Voltage(JROData): """ if channel != None: - data = self.data[channel] + data = self.data[channel,ymin_index:ymax_index] nChannels = 1 else: - data = self.data + data = self.data[:,ymin_index:ymax_index] nChannels = self.nChannels noise = numpy.zeros(nChannels) @@ -407,10 +408,10 @@ class Voltage(JROData): return noise - def getNoise(self, type=1, channel=None): + def getNoise(self, type=1, channel=None,ymin_index=None, ymax_index=None): if type == 1: - noise = self.getNoisebyHildebrand(channel) + noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index) return noise @@ -437,6 +438,8 @@ class Voltage(JROData): class Spectra(JROData): + data_outlier = 0 + def __init__(self): ''' Constructor @@ -475,6 +478,9 @@ class Spectra(JROData): 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles'] + self.max_nIncohInt = 1 + + def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None): """ Determino el nivel de ruido usando el metodo Hildebrand-Sekhon @@ -482,12 +488,28 @@ class Spectra(JROData): Return: noiselevel """ - + # if hasattr(self.nIncohInt, "__len__"): #nIncohInt is a matrix + # + # heis = self.data_spc.shape[2] + # + # noise = numpy.zeros((self.nChannels, heis)) + # for hei in range(heis): + # for channel in range(self.nChannels): + # daux = self.data_spc[channel, xmin_index:xmax_index, hei] + # + # noise[channel,hei] = hildebrand_sekhon(daux, self.nIncohInt[channel,hei]) + # + # else: + # noise = numpy.zeros(self.nChannels) + # for channel in range(self.nChannels): + # daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index] + # + # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt) noise = numpy.zeros(self.nChannels) for channel in range(self.nChannels): - daux = self.data_spc[channel, - xmin_index:xmax_index, ymin_index:ymax_index] - noise[channel] = hildebrand_sekhon(daux, self.nIncohInt) + daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index] + + noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt) return noise @@ -552,6 +574,7 @@ class Spectra(JROData): #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter + return normFactor @property @@ -582,7 +605,7 @@ class Spectra(JROData): def getPower(self): factor = self.normFactor - z = self.data_spc / factor + z = numpy.divide(self.data_spc,factor) z = numpy.where(numpy.isfinite(z), z, numpy.NAN) avg = numpy.average(z, axis=1) diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 2e260c3..fec1ba3 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -58,7 +58,7 @@ 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) + #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) if self.CODE == 'spc_moments': data['moments'] = dataOut.moments @@ -86,9 +86,10 @@ class SpectraPlot(Plot): data = self.data[-1] z = data['spc'] - print(z.shape, x.shape, y.shape) + #print(z.shape, x.shape, y.shape) for n, ax in enumerate(self.axes): - noise = data['noise'][n] + #noise = data['noise'][n] + noise=62 if self.CODE == 'spc_moments': mean = data['moments'][n, 1] if ax.firsttime: @@ -105,15 +106,15 @@ class SpectraPlot(Plot): if self.showprofile: ax.plt_profile = self.pf_axes[n].plot( data['rti'][n], y)[0] - ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y, - color="k", linestyle="dashed", lw=1)[0] + # ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y, + # color="k", linestyle="dashed", lw=1)[0] if self.CODE == 'spc_moments': ax.plt_mean = ax.plot(mean, y, color='k')[0] else: ax.plt.set_array(z[n].T.ravel()) if self.showprofile: ax.plt_profile.set_data(data['rti'][n], y) - ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y) + #ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y) if self.CODE == 'spc_moments': ax.plt_mean.set_data(mean, y) if len(self.azimuthList) > 0 and len(self.elevationList) > 0: @@ -274,7 +275,7 @@ class RTIPlot(Plot): self.x = self.data.times self.y = self.data.yrange - + #print(" x, y: ",self.x, self.y) self.z = self.data[self.CODE] self.z = numpy.array(self.z, dtype=float) self.z = numpy.ma.masked_invalid(self.z) @@ -816,6 +817,7 @@ class NoiselessSpectraPlot(Plot): plot_type = 'pcolor' buffering = False channelList = [] + last_noise = None def setup(self): @@ -842,18 +844,22 @@ class NoiselessSpectraPlot(Plot): self.update_list(dataOut) data = {} meta = {} - n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) - (nch, nff, nh) = dataOut.data_spc.shape - n1 = numpy.repeat(n0,nh, axis=0).reshape((nch,nh)) - noise = numpy.repeat(n1,nff, axis=1).reshape((nch,nff,nh)) - #print(noise.shape, "noise", noise) - spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - noise + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter + n0 = 10*numpy.log10(dataOut.getNoise()/float(norm)) - data['spc'] = spc - data['rti'] = dataOut.getPower() - n1 + 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['noise'] = n0 + data['spc'] = spc - n0 + data['rti'] = dataOut.getPower() - n0 + + #data['noise'] = noise meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) return data, meta @@ -877,7 +883,7 @@ class NoiselessSpectraPlot(Plot): z = data['spc'] for n, ax in enumerate(self.axes): - noise = data['noise'][n] + #noise = data['noise'][n] if ax.firsttime: self.xmax = self.xmax if self.xmax else numpy.nanmax(x) @@ -893,16 +899,15 @@ class NoiselessSpectraPlot(Plot): if self.showprofile: ax.plt_profile = self.pf_axes[n].plot( data['rti'][n], y)[0] - ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y, - color="k", linestyle="dashed", lw=1)[0] + else: ax.plt.set_array(z[n].T.ravel()) if self.showprofile: ax.plt_profile.set_data(data['rti'][n], y) - ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y) - self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise)) + + self.titles.append('CH {}'.format(self.channelList[n])) class NoiselessRTIPlot(Plot): @@ -917,6 +922,7 @@ class NoiselessRTIPlot(Plot): channelList = [] elevationList = [] azimuthList = [] + last_noise = None def setup(self): self.xaxis = 'time' @@ -945,19 +951,21 @@ class NoiselessRTIPlot(Plot): data = {} meta = {} - n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) - (nch, nff, nh) = dataOut.data_spc.shape - #print(nch, nff, nh) - if nch != 1: - aux = [] - for c in self.channelList: - aux.append(n0[c]) - n0 = numpy.asarray(aux) - noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh)) - #print(dataOut.elevationList, dataOut.azimuthList) - #print(dataOut.channelList) - data['noiseless_rti'] = dataOut.getPower() - noise - data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + + 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['noiseless_rti'] = dataOut.getPower() - n0 + + #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor) + #print(noise) return data, meta def plot(self): @@ -986,6 +994,7 @@ class NoiselessRTIPlot(Plot): else: x, y, z = self.fill_gaps(*self.decimate()) 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) @@ -999,9 +1008,6 @@ class NoiselessRTIPlot(Plot): if self.showprofile: ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0] - if "noise" in self.data: - ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y, - color="k", linestyle="dashed", lw=1)[0] else: ax.collections.remove(ax.collections[0]) ax.plt = ax.pcolormesh(x, y, z[n].T, @@ -1011,6 +1017,164 @@ class NoiselessRTIPlot(Plot): ) if self.showprofile: ax.plot_profile.set_data(data['noiseless_rti'][n], self.y) - if "noise" in self.data: - ax.plot_noise.set_data(numpy.repeat( - data['noise'][n], len(self.y)), 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(data['noise'][n], self.y) + + +class OutliersRTIPlot(Plot): + ''' + Plot for data_xxxx object + ''' + + CODE = 'outlier' + colormap = 'cool' + plot_type = 'pcolorbuffer' + + def setup(self): + self.xaxis = 'time' + self.ncols = 1 + self.nrows = self.data.shape('outlier')[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 = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)] + + def update(self, dataOut): + + data = {} + data['outlier'] = dataOut.data_outlier + + 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['outlier'] + + #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['outlier'][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['outlier'][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) diff --git a/schainpy/model/io/jroIO_kamisr.py b/schainpy/model/io/jroIO_kamisr.py index db52c38..d1e722a 100644 --- a/schainpy/model/io/jroIO_kamisr.py +++ b/schainpy/model/io/jroIO_kamisr.py @@ -185,6 +185,7 @@ class AMISRReader(ProcessingUnit): self.nsa = self.nsamplesPulse[0,0] #ngates self.nchannels = len(self.beamCode) self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds + #print("IPPS secs: ",self.ippSeconds) #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 380c4c6..e4a205b 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -20,7 +20,7 @@ from schainpy.model.data.jrodata import hildebrand_sekhon from schainpy.model.data import _noise from schainpy.utils import log - +import matplotlib.pyplot as plt from scipy.optimize import curve_fit class SpectraProc(ProcessingUnit): @@ -125,7 +125,7 @@ class SpectraProc(ProcessingUnit): self.dataOut.flagShiftFFT = False def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False): - + #print("run spc proc") if self.dataIn.type == "Spectra": try: @@ -199,6 +199,7 @@ class SpectraProc(ProcessingUnit): if self.profIndex == nProfiles: self.__updateSpecFromVoltage() + if pairsList == None: self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)] else: @@ -208,6 +209,7 @@ class SpectraProc(ProcessingUnit): self.firstdatatime = None self.profIndex = 0 + else: raise ValueError("The type of input object '%s' is not valid".format( self.dataIn.type)) @@ -529,55 +531,56 @@ class getNoise(Operation): # validacion de velocidades indminPoint = None indmaxPoint = None + if self.dataOut.type == 'Spectra': + if minVel == None and maxVel == None : - if minVel == None and maxVel == None: + freqrange = self.dataOut.getFreqRange(1) - freqrange = self.dataOut.getFreqRange(1) + if minFreq == None: + minFreq = freqrange[0] - if minFreq == None: - minFreq = freqrange[0] + if maxFreq == None: + maxFreq = freqrange[-1] - if maxFreq == None: - maxFreq = freqrange[-1] + if (minFreq < freqrange[0]) or (minFreq > maxFreq): + if self.warnings: + print('minFreq: %.2f is out of the frequency range' % (minFreq)) + print('minFreq is setting to %.2f' % (freqrange[0])) + minFreq = freqrange[0] - if (minFreq < freqrange[0]) or (minFreq > maxFreq): - if self.warnings: - print('minFreq: %.2f is out of the frequency range' % (minFreq)) - print('minFreq is setting to %.2f' % (freqrange[0])) - minFreq = freqrange[0] + if (maxFreq > freqrange[-1]) or (maxFreq < minFreq): + if self.warnings: + print('maxFreq: %.2f is out of the frequency range' % (maxFreq)) + print('maxFreq is setting to %.2f' % (freqrange[-1])) + maxFreq = freqrange[-1] - if (maxFreq > freqrange[-1]) or (maxFreq < minFreq): - if self.warnings: - print('maxFreq: %.2f is out of the frequency range' % (maxFreq)) - print('maxFreq is setting to %.2f' % (freqrange[-1])) - maxFreq = freqrange[-1] + indminPoint = numpy.where(freqrange >= minFreq) + indmaxPoint = numpy.where(freqrange <= maxFreq) - indminPoint = numpy.where(freqrange >= minFreq) - indmaxPoint = numpy.where(freqrange <= maxFreq) + else: - else: - velrange = self.dataOut.getVelRange(1) + velrange = self.dataOut.getVelRange(1) - if minVel == None: - minVel = velrange[0] + if minVel == None: + minVel = velrange[0] - if maxVel == None: - maxVel = velrange[-1] + if maxVel == None: + maxVel = velrange[-1] - if (minVel < velrange[0]) or (minVel > maxVel): - if self.warnings: - print('minVel: %.2f is out of the velocity range' % (minVel)) - print('minVel is setting to %.2f' % (velrange[0])) - minVel = velrange[0] + if (minVel < velrange[0]) or (minVel > maxVel): + if self.warnings: + print('minVel: %.2f is out of the velocity range' % (minVel)) + print('minVel is setting to %.2f' % (velrange[0])) + minVel = velrange[0] - if (maxVel > velrange[-1]) or (maxVel < minVel): - if self.warnings: - print('maxVel: %.2f is out of the velocity range' % (maxVel)) - print('maxVel is setting to %.2f' % (velrange[-1])) - maxVel = velrange[-1] + if (maxVel > velrange[-1]) or (maxVel < minVel): + if self.warnings: + print('maxVel: %.2f is out of the velocity range' % (maxVel)) + print('maxVel is setting to %.2f' % (velrange[-1])) + maxVel = velrange[-1] - indminPoint = numpy.where(velrange >= minVel) - indmaxPoint = numpy.where(velrange <= maxVel) + indminPoint = numpy.where(velrange >= minVel) + indmaxPoint = numpy.where(velrange <= maxVel) # seleccion de indices para rango @@ -606,23 +609,30 @@ class getNoise(Operation): maxIndex = self.dataOut.nHeights - 1 #############################################################3 # seleccion de indices para velocidades + if self.dataOut.type == 'Spectra': + try: + minIndexFFT = indminPoint[0][0] + except: + minIndexFFT = 0 - try: - minIndexFFT = indminPoint[0][0] - except: - minIndexFFT = 0 + try: + maxIndexFFT = indmaxPoint[0][-1] + except: + maxIndexFFT = len( self.dataOut.getFreqRange(1)) - try: - maxIndexFFT = indmaxPoint[0][-1] - except: - maxIndexFFT = len( self.dataOut.getFreqRange(1)) - #print(minIndex, maxIndex,minIndexVel, maxIndexVel) self.dataOut.noise_estimation = None - noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex) - + noise = None + if self.dataOut.type == 'Voltage': + noise = self.dataOut.getNoise(ymin_index=minIndex, ymax_index=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) + else: + noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex) self.dataOut.noise_estimation = noise.copy() # dataOut.noise #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64)) + return self.dataOut @@ -705,7 +715,7 @@ class CleanRayleigh(Operation): def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5): - + #print("runing cleanRayleigh") if not self.isConfig : self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv) @@ -736,7 +746,10 @@ class CleanRayleigh(Operation): if numpy.any(jspc) : #print( jspc.shape, jcspc.shape) jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights)) - jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights)) + try: + jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights)) + except: + print("no cspc") self.__dataReady = False #print( jspc.shape, jcspc.shape) dataOut.flagNoData = False @@ -748,8 +761,11 @@ class CleanRayleigh(Operation): #print( len(self.buffer)) if numpy.any(self.buffer): self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0) - self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0) - self.buffer3 += dataOut.data_dc + try: + self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0) + self.buffer3 += dataOut.data_dc + except: + pass else: self.buffer = dataOut.data_spc self.buffer2 = dataOut.data_cspc @@ -762,7 +778,9 @@ 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 @@ -776,6 +794,7 @@ class CleanRayleigh(Operation): dataOut.data_dc = self.buffer3 dataOut.nIncohInt *= self.nIntProfiles + dataOut.max_nIncohInt = self.nIntProfiles dataOut.utctime = self.currentTime #tiempo promediado #print("Time: ",time.localtime(dataOut.utctime)) # dataOut.data_spc = sat_spectra @@ -787,7 +806,7 @@ class CleanRayleigh(Operation): return dataOut def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv): - #print("OP cleanRayleigh") + print("OP cleanRayleigh") #import matplotlib.pyplot as plt #for k in range(149): #channelsProcssd = [] @@ -803,6 +822,8 @@ class CleanRayleigh(Operation): ###ONLY FOR TEST: raxs = math.ceil(math.sqrt(self.nPairs)) + if raxs == 0: + raxs = 1 caxs = math.ceil(self.nPairs/raxs) if self.nPairs <4: raxs = 2 @@ -1100,12 +1121,13 @@ class IntegrationFaradaySpectra(Operation): __dataReady = False __timeInterval = None - + n_ints = None #matriz de numero de integracions (CH,HEI) n = None minHei_ind = None maxHei_ind = None navg = 1.0 factor = 0.0 + dataoutliers = None # (CHANNELS, HEIGHTS) def __init__(self): @@ -1138,6 +1160,7 @@ class IntegrationFaradaySpectra(Operation): self.navg = avg #self.ByLags = dataOut.ByLags ###REDEFINIR self.ByLags = False + self.maxProfilesInt = 1 if DPL != None: self.DPL=DPL @@ -1251,6 +1274,8 @@ class IntegrationFaradaySpectra(Operation): freq_dc = int(self.__buffer_spc.shape[2] / 2) #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights) + self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers + for k in range(self.minHei_ind,self.maxHei_ind): try: buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k]) @@ -1266,22 +1291,31 @@ class IntegrationFaradaySpectra(Operation): #sortIDs=[] outliers_IDs=[] - for j in range(self.nProfiles): + for j in range(self.nProfiles): #frecuencias en el tiempo # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0 # continue # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1 # continue buffer=buffer1[:,j] sortdata = numpy.sort(buffer, axis=None) + sortID=buffer.argsort() index = _noise.hildebrand_sekhon2(sortdata,self.navg) #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor) + # fig,ax = plt.subplots() + # ax.set_title(str(k)+" "+str(j)) + # x=range(len(sortdata)) + # ax.scatter(x,sortdata) + # ax.axvline(index) + # plt.show() + indexes.append(index) #sortIDs.append(sortID) outliers_IDs=numpy.append(outliers_IDs,sortID[index:]) + #print("Outliers: ",outliers_IDs) outliers_IDs=numpy.array(outliers_IDs) outliers_IDs=outliers_IDs.ravel() outliers_IDs=numpy.unique(outliers_IDs) @@ -1289,31 +1323,45 @@ class IntegrationFaradaySpectra(Operation): indexes=numpy.array(indexes) indexmin=numpy.min(indexes) + + #print(indexmin,buffer1.shape[0], k) + + # fig,ax = plt.subplots() + # ax.plot(sortdata) + # ax2 = ax.twinx() + # x=range(len(indexes)) + # #plt.scatter(x,indexes) + # ax2.scatter(x,indexes) + # plt.show() + if indexmin != buffer1.shape[0]: if self.nChannels > 1: cspc_outliers_exist= True lt=outliers_IDs - avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0) + #avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0) for p in list(outliers_IDs): - buffer1[p,:]=avg + #buffer1[p,:]=avg + buffer1[p,:] = numpy.NaN + + self.dataOutliers[i,k] = len(outliers_IDs) self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1) - ###cspc IDs - #indexmin_cspc+=indexmin_cspc + outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs) - #if not breakFlag: + + outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64')) if cspc_outliers_exist : - #sortdata=numpy.sort(buffer_cspc,axis=0) - #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0) + lt=outliers_IDs_cspc - avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0) + #avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0) for p in list(outliers_IDs_cspc): - buffer_cspc[p,:]=avg + #buffer_cspc[p,:]=avg + buffer_cspc[p,:] = numpy.NaN try: self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc) @@ -1325,39 +1373,36 @@ class IntegrationFaradaySpectra(Operation): - + nOutliers = len(outliers_IDs) + #print("Outliers n: ",self.dataOutliers,nOutliers) buffer=None bufferH=None buffer1=None buffer_cspc=None - #print("cpsc",self.__buffer_cspc[:,0,0,0,0]) - #exit() buffer=None - #print(self.__buffer_spc[:,1,3,20,0]) - #print(self.__buffer_spc[:,1,5,37,0]) - data_spc = numpy.sum(self.__buffer_spc,axis=0) + + #data_spc = numpy.sum(self.__buffer_spc,axis=0) + data_spc = numpy.nansum(self.__buffer_spc,axis=0) try: - data_cspc = numpy.sum(self.__buffer_cspc,axis=0) + #data_cspc = numpy.sum(self.__buffer_cspc,axis=0) + data_cspc = numpy.nansum(self.__buffer_cspc,axis=0) except: #print("No cpsc",e) pass - #print(numpy.shape(data_spc)) - #data_spc[1,4,20,0]=numpy.nan - - #data_cspc = self.__buffer_cspc - #print("pushData pre Done") data_dc = self.__buffer_dc - n = self.__profIndex + #(CH, HEIGH) + self.maxProfilesInt = self.__profIndex + n = self.__profIndex - self.dataOutliers self.__buffer_spc = [] self.__buffer_cspc = [] self.__buffer_dc = 0 self.__profIndex = 0 - #print("pushData Done") + return data_spc, data_cspc, data_dc, n def byProfiles(self, *args): @@ -1372,7 +1417,7 @@ class IntegrationFaradaySpectra(Operation): if self.__profIndex == self.n: avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData() - self.n = n + self.n_ints = n self.__dataReady = True return avgdata_spc, avgdata_cspc, avgdata_dc @@ -1388,7 +1433,7 @@ class IntegrationFaradaySpectra(Operation): if (datatime - self.__initime) >= self.__integrationtime: avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData() - self.n = n + self.n_ints = n self.__dataReady = True return avgdata_spc, avgdata_cspc, avgdata_dc @@ -1450,7 +1495,7 @@ class IntegrationFaradaySpectra(Operation): self.dataOut.data_spc = numpy.squeeze(avgdata_spc) self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc) self.dataOut.data_dc = avgdata_dc - + self.dataOut.data_outlier = self.dataOutliers else: self.dataOut.dataLag_spc = avgdata_spc @@ -1462,10 +1507,12 @@ class IntegrationFaradaySpectra(Operation): self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot] - self.dataOut.nIncohInt *= self.n + self.dataOut.nIncohInt *= self.n_ints + self.dataOut.max_nIncohInt = self.maxProfilesInt + #print(self.dataOut.max_nIncohInt) self.dataOut.utctime = avgdatatime self.dataOut.flagNoData = False - + #print("Faraday Integration DONE...") return self.dataOut class removeInterference(Operation): @@ -1705,7 +1752,8 @@ class IncohInt(Operation): __dataReady = False __timeInterval = None - + incohInt = 0 + nOutliers = 0 n = None def __init__(self): @@ -1734,7 +1782,8 @@ class IncohInt(Operation): self.__profIndex = 0 self.__dataReady = False self.__byTime = False - + self.incohInt = 0 + self.nOutliers = 0 if n is None and timeInterval is None: raise ValueError("n or timeInterval should be specified ...") @@ -1788,7 +1837,7 @@ class IncohInt(Operation): self.__buffer_spc = 0 self.__buffer_cspc = 0 self.__buffer_dc = 0 - self.__profIndex = 0 + return data_spc, data_cspc, data_dc, n @@ -1845,6 +1894,9 @@ class IncohInt(Operation): if n == 1: return dataOut + if dataOut.flagNoData == True: + return dataOut + dataOut.flagNoData = True if not self.isConfig: @@ -1855,15 +1907,21 @@ class IncohInt(Operation): dataOut.data_spc, dataOut.data_cspc, dataOut.data_dc) - + self.incohInt += dataOut.nIncohInt + self.nOutliers += dataOut.data_outlier if self.__dataReady: dataOut.data_spc = avgdata_spc dataOut.data_cspc = avgdata_cspc dataOut.data_dc = avgdata_dc - dataOut.nIncohInt *= self.n + dataOut.nIncohInt = self.incohInt + dataOut.data_outlier = self.nOutliers dataOut.utctime = avgdatatime dataOut.flagNoData = False + dataOut.max_nIncohInt = self.__profIndex + self.incohInt = 0 + self.nOutliers = 0 + self.__profIndex = 0 return dataOut diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 0e821a4..2d59d67 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -8,7 +8,8 @@ from schainpy.model.io.utils import getHei_index from time import time import datetime import numpy -import copy +#import copy +from schainpy.model.data import _noise class VoltageProc(ProcessingUnit): @@ -1759,6 +1760,8 @@ class SSheightProfiles2(Operation): profileShape = None sshProfiles = None profileIndex = None + deltaHeight = None + init_range = None def __init__(self, **kwargs): @@ -1780,7 +1783,8 @@ class SSheightProfiles2(Operation): if residue != 0: print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue)) - deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] + 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 @@ -1800,21 +1804,25 @@ class SSheightProfiles2(Operation): if repeat is not None: code_block = numpy.repeat(code_block, repeats=repeat, axis=1) - #print("buff, data, :",self.buffer.shape, data.shape) + if data.ndim == 2: + data = data.reshape(1,1,self.__nHeis ) + #print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape) for i in range(int(self.new_nHeights)): #nuevas alturas if code is not None: - self.buffer[:,i] = data[:,i*self.step:i*self.step + self.nsamples]*code_block + self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]*code_block else: - self.buffer[:,i] = data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:] + self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:] for j in range(self.__nChannels): #en los cananles - self.sshProfiles[j] = numpy.transpose(self.buffer[j]) - + self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:]) + #print("new profs Done") def run(self, dataOut, step, nsamples, code = None, repeat = None): + if dataOut.flagNoData == True: + return dataOut dataOut.flagNoData = True #print("init data shape:", dataOut.data.shape) #print("ch: {} prof: {} hs: {}".format(int(dataOut.nChannels), @@ -1838,24 +1846,26 @@ class SSheightProfiles2(Operation): #print("dataOut nProfiles:", dataOut.nProfiles) for profile in range(nprof): if dataOut.flagDataAsBlock: + #print("read blocks") self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat) else: - self.getNewProfiles(dataOut.data, code=code, repeat=repeat) + #print("read profiles") + self.getNewProfiles(dataOut.data, code=code, repeat=repeat) #only one channe if profile == 0: dataBlock = self.sshProfiles.copy() else: #by blocks dataBlock = numpy.concatenate((dataBlock,self.sshProfiles), axis=1) #profile axis - print("by blocks: ",dataBlock.shape, self.sshProfiles.shape) + #print("by blocks: ",dataBlock.shape, self.sshProfiles.shape) profileIndex = self.nsamples - deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] - ippSeconds = (deltaHeight*1.0e-6)/(0.15) + #deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] + ippSeconds = (self.deltaHeight*1.0e-6)/(0.15) dataOut.data = dataBlock - dataOut.heightList = numpy.arange(self.new_nHeights) *self.step*deltaHeight + dataOut.heightList[0] - #print("show me: ",dataOut.nHeights, self.new_nHeights) - #dataOut.nHeights = int(self.new_nHeights) + #print("show me: ",self.step,self.deltaHeight, dataOut.heightList, self.new_nHeights) + dataOut.heightList = numpy.arange(int(self.new_nHeights)) *self.step*self.deltaHeight + self.init_range + dataOut.ippSeconds = ippSeconds dataOut.step = self.step dataOut.flagNoData = False @@ -1877,17 +1887,17 @@ class SSheightProfiles2(Operation): #import skimage.color #import skimage.io -import matplotlib.pyplot as plt +#import matplotlib.pyplot as plt -class clean2DArray(Operation): +class removeProfileByFaradayHS(Operation): ''' ''' isConfig = False n = None - __timeInterval = None + __profIndex = 0 - __byTime = False + __dataReady = False __buffer_data = [] __buffer_times = [] @@ -1901,6 +1911,7 @@ class clean2DArray(Operation): end_prof = 0 n_profiles = 0 first_utcBlock = None + outliers_IDs_list = [] __dh = 0 def __init__(self, **kwargs): @@ -1908,42 +1919,122 @@ class clean2DArray(Operation): Operation.__init__(self, **kwargs) self.isConfig = False - def setup(self,dataOut, n=None , timeInterval=None): + def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=3, minHei=None, maxHei=None): if n == None and timeInterval == None: raise ValueError("nprofiles or timeInterval should be specified ...") if n != None: self.n = n - self.__byTime = False - else: - self.__timeInterval = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line - self.n = 9999 - self.__byTime = True + self.navg = navg + self.profileMargin = profileMargin + self.thHistOutlier = thHistOutlier self.__profIndex = 0 self.buffer = None self.lenProfileOut = 1 - + self.n_prof_released = 0 + self.heightList = dataOut.heightList self.init_prof = 0 self.end_prof = 0 self.n_profiles = 0 self.first_utcBlock = None self.__dh = dataOut.heightList[1] - dataOut.heightList[0] + minHei = minHei + maxHei = maxHei + if minHei==None : + minHei = dataOut.heightList[0] + if maxHei==None : + maxHei = dataOut.heightList[-1] + self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList) + + + def filterSatsProfiles(self): + data = self.__buffer_data + #print(data.shape) + nChannels, profiles, heights = data.shape + indexes=[] + outliers_IDs=[] + for c in range(nChannels): + for h in range(self.minHei_idx, self.maxHei_idx): + power = data[c,:,h] * numpy.conjugate(data[c,:,h]) + power = power.real + #power = (numpy.abs(data[c,:,h].real)) + sortdata = numpy.sort(power, axis=None) + sortID=power.argsort() + index = _noise.hildebrand_sekhon2(sortdata,self.navg) #0.75-> buen valor + + indexes.append(index) + outliers_IDs=numpy.append(outliers_IDs,sortID[index:]) + # print(outliers_IDs) + # fig,ax = plt.subplots() + # #ax.set_title(str(k)+" "+str(j)) + # x=range(len(sortdata)) + # ax.scatter(x,sortdata) + # ax.axvline(index) + # plt.grid() + # plt.show() + + outliers_IDs = outliers_IDs.astype(numpy.dtype('int64')) + outliers_IDs = numpy.unique(outliers_IDs) + outs_lines = numpy.sort(outliers_IDs) + # #print("outliers Ids: ", outs_lines, outs_lines.shape) + #hist, bin_edges = numpy.histogram(outs_lines, bins=10, density=True) + + + #Agrupando el histograma de outliers, + my_bins = numpy.linspace(0,9600, 96, endpoint=False) + + + hist, bins = numpy.histogram(outs_lines,bins=my_bins) + hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier + #print(hist_outliers_indexes[0]) + bins_outliers_indexes = [int(i) for i in bins[hist_outliers_indexes]] # + #print(bins_outliers_indexes) + outlier_loc_index = [] + + #outlier_loc_index = [k for k in range(bins_outliers_indexes[n]-50,bins_outliers_indexes[n+1]+50) for n in range(len(bins_outliers_indexes)-1) ] + for n in range(len(bins_outliers_indexes)-1): + #outlier_loc_index = [k for k in range(bins_outliers_indexes[n]-50,bins_outliers_indexes[n+1]+50)] + for k in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin): + outlier_loc_index.append(k) + + outlier_loc_index = numpy.asarray(outlier_loc_index) + #print(numpy.unique(outlier_loc_index)) + + + + # x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList) + # fig, ax = plt.subplots(1,2,figsize=(8, 6)) + # + # dat = data[0,:,:].real + # m = numpy.nanmean(dat) + # o = numpy.nanstd(dat) + # #print(m, o, x.shape, y.shape) + # c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o)) + # ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w') + # fig.colorbar(c) + # ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r') + # ax[1].hist(outs_lines,bins=my_bins) + # plt.show() + + + self.outliers_IDs_list = numpy.unique(outlier_loc_index) + return data def cleanOutliersByBlock(self): #print(self.__buffer_data[0].shape) data = self.__buffer_data#.copy() #print("cleaning shape inpt: ",data.shape) - + ''' self.__buffer_data = [] spectrum = numpy.fft.fft2(data, axes=(0,2)) #print("spc : ",spectrum.shape) (nch,nsamples, nh) = spectrum.shape - - #nch = - + data2 = None + #print(data.shape) + spectrum2 = spectrum.copy() for ch in range(nch): dh = self.__dh dt1 = (dh*1.0e-6)/(0.15) @@ -1951,57 +2042,48 @@ class clean2DArray(Operation): freqv = numpy.fft.fftfreq(nh, d=dt1) freqh = numpy.fft.fftfreq(self.n, d=dt2) #print("spc loop: ") - x, y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv)) - z = numpy.abs(spectrum[ch,:,:]) - - - dat = numpy.log10(z.T) - m = numpy.mean(dat) - o = numpy.std(dat) - fig, ax = plt.subplots(figsize=(8, 6)) - c = ax.pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o)) - #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o)) - date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f') - #strftime('%Y-%m-%d %H:%M:%S') - ax.set_title('Spectrum magnitude '+date_time) - fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time) - fig.colorbar(c) - plt.show() + x, y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv)) + z = numpy.abs(spectrum[ch,:,:]) + # Find all peaks higher than the 98th percentile + peaks = z < numpy.percentile(z, 98) + #print(peaks) + # Set those peak coefficients to zero + spectrum2 = spectrum2 * peaks.astype(int) + data2 = numpy.fft.ifft2(spectrum2) + dat = numpy.log10(z.T) + dat2 = numpy.log10(spectrum2.T) + + # m = numpy.mean(dat) + # o = numpy.std(dat) + # fig, ax = plt.subplots(2,1,figsize=(8, 6)) + # + # c = ax[0].pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o)) + # #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o)) + # date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f') + # #strftime('%Y-%m-%d %H:%M:%S') + # ax[0].set_title('Spectrum magnitude '+date_time) + # fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time) + # + # + # c = ax[1].pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o)) + # fig.colorbar(c) + # plt.show() + + #print(data2.shape) + + data = data2 #cleanBlock = numpy.fft.ifft2(spectrum, axes=(0,2)).reshape() + ''' + #print("cleanOutliersByBlock Done") - print("cleanOutliersByBlock Done") - return data - - def byTime(self, data, datatime): - - self.__dataReady = False - - self.fillBuffer(data, datatime) - dataBlock = None - - if (datatime - self.__initime) >= self.__timeInterval: - dataBlock = self.cleanOutliersByBlock() - self.__dataReady = True - self.n = self.__profIndex - return dataBlock - - def byProfiles(self, data, datatime): - self.__dataReady = False - - self.fillBuffer(data, datatime) - dataBlock = None + return self.filterSatsProfiles() - if self.__profIndex == self.n: - #print("apnd : ",data) - dataBlock = self.cleanOutliersByBlock() - self.__dataReady = True - return dataBlock def fillBuffer(self, data, datatime): @@ -2011,54 +2093,60 @@ class clean2DArray(Operation): else: self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles self.__profIndex += 1 - self.__buffer_times.append(datatime) + #self.__buffer_times.append(datatime) def getData(self, data, datatime=None): if self.__profIndex == 0: self.__initime = datatime - if self.__byTime: - dataBlock = self.byTime(data, datatime) - else: - dataBlock = self.byProfiles(data, datatime) + self.__dataReady = False + self.fillBuffer(data, datatime) + dataBlock = None + + if self.__profIndex == self.n: + #print("apnd : ",data) + #dataBlock = self.cleanOutliersByBlock() + dataBlock = self.filterSatsProfiles() + self.__dataReady = True + + return dataBlock if dataBlock is None: - return None + return None, None return dataBlock - def releaseBlock(self, dataOut): + def releaseBlock(self): if self.n % self.lenProfileOut != 0: raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles)) return None - dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt - dataOut.profileIndex = self.lenProfileOut - dataOut.utctime = self.first_utcBlock + self.init_prof*dataOut.ippSeconds + data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt + self.init_prof = self.end_prof self.end_prof += self.lenProfileOut - print("data release shape: ",dataOut.data.shape, self.end_prof) + #print("data release shape: ",dataOut.data.shape, self.end_prof) + self.n_prof_released += 1 - if self.end_prof >= (self.n +self.lenProfileOut): - self.init_prof = 0 - self.__profIndex = 0 - self.buffer = None - dataOut.buffer_empty = True - return dataOut - def run(self, dataOut, n=None, timeInterval=None, nProfilesOut=1): + #print("f_no_data ", dataOut.flagNoData) + 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 if not self.isConfig: - self.setup(dataOut,n=n, timeInterval=timeInterval) + #print("init p idx: ", dataOut.profileIndex ) + self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin, + thHistOutlier=th_hist_outlier,minHei=minHei, maxHei=maxHei) self.isConfig = True dataBlock = None @@ -2066,6 +2154,7 @@ class clean2DArray(Operation): if not dataOut.buffer_empty: #hay datos acumulados if self.init_prof == 0: + self.n_prof_released = 0 self.lenProfileOut = nProfilesOut dataOut.flagNoData = False #print("tp 2 ",dataOut.data.shape) @@ -2074,13 +2163,45 @@ class clean2DArray(Operation): #print("rel: ",self.buffer[:,-1,:]) self.init_prof = 0 self.end_prof = self.lenProfileOut - self.first_utcBlock = dataOut.utctime + dataOut.nProfiles = self.lenProfileOut - dataOut.error = False + if nProfilesOut == 1: + dataOut.flagDataAsBlock = False + else: dataOut.flagDataAsBlock = True #print("prof: ",self.init_prof) dataOut.flagNoData = False - return self.releaseBlock(dataOut) + 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.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)) + return None + + dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt + + self.init_prof = self.end_prof + self.end_prof += self.lenProfileOut + #print("data release shape: ",dataOut.data.shape, self.end_prof) + self.n_prof_released += 1 + + if self.end_prof >= (self.n +self.lenProfileOut): + + self.init_prof = 0 + self.__profIndex = 0 + self.buffer = None + dataOut.buffer_empty = True + self.outliers_IDs_list = [] + self.n_prof_released = 0 + dataOut.flagNoData = False #enviar ultimo aunque sea outlier :( + #print("cleaning...", dataOut.buffer_empty) + dataOut.profileIndex = 0 #self.lenProfileOut + #################################################################### + return dataOut #print("tp 223 ",dataOut.data.shape) @@ -2098,13 +2219,16 @@ class clean2DArray(Operation): sys.exit() if self.__dataReady: + #print("omitting: ", len(self.outliers_IDs_list)) self.__count_exec = 0 #dataOut.data = - self.buffer = numpy.flip(dataBlock, axis=1) - + #self.buffer = numpy.flip(dataBlock, axis=1) + self.buffer = dataBlock + self.first_utcBlock = self.__initime dataOut.utctime = self.__initime dataOut.nProfiles = self.__profIndex #dataOut.flagNoData = False + self.init_prof = 0 self.__profIndex = 0 self.__initime = None dataBlock = None @@ -2112,16 +2236,15 @@ class clean2DArray(Operation): dataOut.error = False dataOut.useInputBuffer = True dataOut.buffer_empty = False - print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels), - int(dataOut.nProfiles),int(dataOut.nHeights))) - #return None + #print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights))) + #print(self.__count_exec) return dataOut -class CleanProfileSats(Operation): +class RemoveProfileSats(Operation): ''' Omite los perfiles contaminados con señal de satelites, In: minHei = min_sat_range @@ -2144,6 +2267,7 @@ class CleanProfileSats(Operation): byRanges = False min_sats = None max_sats = None + noise = 0 def __init__(self, **kwargs): @@ -2177,16 +2301,19 @@ class CleanProfileSats(Operation): def compareRanges(self,data, minHei,maxHei): - ref = data[0,self.min_ref:self.max_ref] * numpy.conjugate(data[0,self.min_ref:self.max_ref]) - p_ref = 10*numpy.log10(ref.real) - m_ref = numpy.mean(p_ref) + + # ref = data[0,self.min_ref:self.max_ref] * numpy.conjugate(data[0,self.min_ref:self.max_ref]) + # p_ref = 10*numpy.log10(ref.real) + # m_ref = numpy.mean(p_ref) + + m_ref = self.noise sats = data[0,minHei:maxHei] * numpy.conjugate(data[0,minHei:maxHei]) p_sats = 10*numpy.log10(sats.real) m_sats = numpy.mean(p_sats) - if m_sats > (m_ref + self.th) and (m_sats > self.thdB): - #print("msats: ",m_sats," mRef: ", m_ref, (m_sats - m_ref)) + if m_sats > (m_ref + self.th): #and (m_sats > self.thdB): + #print("msats: ",m_sats," \tmRef: ", m_ref, "\t",(m_sats - m_ref)) #print("Removing profiles...") return False @@ -2218,12 +2345,12 @@ class CleanProfileSats(Operation): if not self.isConfig: self.setup(dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList) self.isConfig = True - + #print(self.min_sats,self.max_sats) if dataOut.flagDataAsBlock: raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False") else: - + self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref)) if not self.isProfileClean(dataOut.data): return dataOut #dataOut.data = numpy.full((dataOut.nChannels,dataOut.nHeights),numpy.NAN)