diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 78f3838..b30bb08 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -60,7 +60,7 @@ class SpectraPlot(Plot): data['rti'] = dataOut.getPower() norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter noise = 10*numpy.log10(dataOut.getNoise()/float(norm)) - data['noise'] = noise[0] + data['noise'] = noise meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) if self.CODE == 'spc_moments': @@ -91,7 +91,7 @@ class SpectraPlot(Plot): z = data['spc'] #print(z.shape, x.shape, y.shape) for n, ax in enumerate(self.axes): - noise = self.data['noise'][n] + noise = self.data['noise'][n][0] #print(noise) if self.CODE == 'spc_moments': mean = data['moments'][n, 1] @@ -275,7 +275,7 @@ class RTIPlot(Plot): 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): @@ -533,9 +533,10 @@ class SpectraCutPlot(Plot): 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 + spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) + noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights) - data['spc'] = spc + data['spc'] = spc - noise meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) return data, meta @@ -827,7 +828,7 @@ class NoiselessSpectraPlot(Plot): ''' CODE = 'noiseless_spc' - colormap = 'nipy_spectral' + colormap = 'jet' plot_type = 'pcolor' buffering = False channelList = [] @@ -866,8 +867,12 @@ class NoiselessSpectraPlot(Plot): spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - data['spc'] = spc - n0 - data['rti'] = dataOut.getPower() - n0 + noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights) + data['rti'] = dataOut.getPower() - noise + + noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights) + data['spc'] = spc - noise + # data['noise'] = noise meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) @@ -963,11 +968,13 @@ class NoiselessRTIPlot(Plot): norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter + #print("Norm: ", norm) + #print(dataOut.nProfiles , dataOut.max_nIncohInt ,dataOut.nCohInt, dataOut.windowOfFilter) n0 = 10*numpy.log10(dataOut.getNoise()/float(norm)) - data['noise'] = n0[0] - - data['noiseless_rti'] = dataOut.getPower() - n0 + data['noise'] = n0 + noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights) + data['noiseless_rti'] = dataOut.getPower() - noise return data, meta diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index b63112d..14c0346 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -128,7 +128,7 @@ class ParametersProc(ProcessingUnit): #---------------------- Spectra Data --------------------------- if self.dataIn.type == "Spectra": - + 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 diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 9f0c31d..d884ad6 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -22,6 +22,7 @@ from schainpy.model.data import _noise from schainpy.utils import log import matplotlib.pyplot as plt #from scipy.optimize import curve_fit +from schainpy.model.io.utils import getHei_index class SpectraProc(ProcessingUnit): @@ -73,6 +74,7 @@ class SpectraProc(ProcessingUnit): def __getFft(self): + # print("fft donw") """ Convierte valores de Voltaje a Spectra @@ -124,7 +126,7 @@ class SpectraProc(ProcessingUnit): self.dataOut.blockSize = blocksize self.dataOut.flagShiftFFT = False - def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False): + def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, zeroPad=False): #print("run spc proc") try: type = self.dataIn.type.decode("utf-8") @@ -167,15 +169,21 @@ class SpectraProc(ProcessingUnit): self.dataOut.nFFTPoints = nFFTPoints #print(" volts ch,prof, h: ", self.dataIn.data.shape) if self.buffer is None: - self.buffer = numpy.zeros((self.dataIn.nChannels, + if not zeroPad: + self.buffer = numpy.zeros((self.dataIn.nChannels, nProfiles, self.dataIn.nHeights), dtype='complex') + else: + self.buffer = numpy.zeros((self.dataIn.nChannels, + nFFTPoints, + self.dataIn.nHeights), + dtype='complex') if self.dataIn.flagDataAsBlock: nVoltProfiles = self.dataIn.data.shape[1] - if nVoltProfiles == nProfiles: + if nVoltProfiles == nProfiles or zeroPad: self.buffer = self.dataIn.data.copy() self.profIndex = nVoltProfiles @@ -201,7 +209,7 @@ class SpectraProc(ProcessingUnit): if self.firstdatatime == None: self.firstdatatime = self.dataIn.utctime - if self.profIndex == nProfiles: + if self.profIndex == nProfiles or zeroPad: self.__updateSpecFromVoltage() @@ -682,6 +690,7 @@ class getNoiseB(Operation): #print(self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.nIncohInt) noise = numpy.zeros( self.dataOut.nChannels) norm = 1 + for channel in range( self.dataOut.nChannels): if not hasattr(self.dataOut.nIncohInt,'__len__'): norm = 1 @@ -691,10 +700,11 @@ class getNoiseB(Operation): 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.getNoiseByMean(daux)/self.offset #print(daux.shape, daux) - noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset - + #noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset + sortdata = numpy.sort(daux, axis=None) + noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt)/self.offset # data = numpy.mean(daux,axis=1) # sortdata = numpy.sort(data, axis=None) # noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt)/self.offset @@ -743,7 +753,7 @@ class getNoiseB(Operation): sortdata = numpy.sort(data, axis=None) lenOfData = len(sortdata) - nums_min = lenOfData*0.05 + nums_min = lenOfData*0.2 if nums_min <= 5: @@ -761,7 +771,7 @@ class getNoiseB(Operation): sumq += sortdata[j]**2 #sumq -= sump**2 if j > nums_min: - rtest = float(j)/(j-1) + 1.0/0.1 + rtest = float(j)/(j-1) + 1.0/navg #if ((sumq*j) > (sump**2)): if ((sumq*j) > (rtest*sump**2)): j = j - 1 @@ -1648,8 +1658,74 @@ class IntegrationFaradaySpectra(Operation): #print(self.dataOut.flagNoData) return self.dataOut + + class removeInterference(Operation): + def removeInterference3(self, min_hei = None, max_hei = None): + + jspectra = self.dataOut.data_spc + #jcspectra = self.dataOut.data_cspc + jnoise = self.dataOut.getNoise() + num_incoh = self.dataOut.max_nIncohInt + #print(jspectra.shape) + num_channel, num_prof, num_hei = jspectra.shape + minHei = min_hei + maxHei = max_hei + ######################################################################## + if minHei == None or (minHei < self.dataOut.heightList[0]): + minHei = self.dataOut.heightList[0] + + if maxHei == None or (maxHei > self.dataOut.heightList[-1]): + maxHei = self.dataOut.heightList[-1] + minIndex = 0 + maxIndex = 0 + heights = self.dataOut.heightList + + inda = numpy.where(heights >= minHei) + indb = numpy.where(heights <= maxHei) + + try: + minIndex = inda[0][0] + except: + minIndex = 0 + try: + maxIndex = indb[0][-1] + except: + maxIndex = len(heights) + + if (minIndex < 0) or (minIndex > maxIndex): + raise ValueError("some value in (%d,%d) is not valid" % ( + minIndex, maxIndex)) + if (maxIndex >= self.dataOut.nHeights): + maxIndex = self.dataOut.nHeights - 1 + + ######################################################################## + + + #dataOut.max_nIncohInt * dataOut.nCohInt + norm = self.dataOut.nIncohInt /self.dataOut.max_nIncohInt + #print(norm.shape) + # Subrutina de Remocion de la Interferencia + for ich in range(num_channel): + # Se ordena los espectros segun su potencia (menor a mayor) + #power = jspectra[ich, mask_prof, :] + interf = jspectra[ich, :, minIndex:maxIndex] + #print(interf.shape) + inttef = interf.mean(axis=1) + + for hei in range(num_hei): + temp = jspectra[ich,:, hei] + temp -= inttef + temp += jnoise[ich]*norm[ich,hei] + jspectra[ich,:, hei] = temp + + # Guardar Resultados + self.dataOut.data_spc = jspectra + #self.dataOut.data_cspc = jcspectra + + return 1 + def removeInterference2(self): cspc = self.dataOut.data_cspc @@ -1687,7 +1763,7 @@ class removeInterference(Operation): # hei_interf if hei_interf is None: - count_hei = int(num_hei / 2) + count_hei = int(num_hei / 2) # a half of total ranges hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei hei_interf = numpy.asarray(hei_interf)[0] #print(hei_interf) @@ -1720,7 +1796,7 @@ class removeInterference(Operation): power = power[:, hei_interf] power = power.sum(axis=0) psort = power.ravel().argsort() - + print(hei_interf[psort[list(range(offhei_interf, nhei_interf + offhei_interf))]]) # Se estima la interferencia promedio en los Espectros de Potencia empleando junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range( offhei_interf, nhei_interf + offhei_interf))]]] @@ -1730,7 +1806,7 @@ class removeInterference(Operation): tmp_noise = jnoise[ich] junkspc_interf = junkspc_interf - tmp_noise #junkspc_interf[:,comp_mask_prof] = 0 - + print(junkspc_interf.shape) jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf jspc_interf = jspc_interf.transpose() # Calculando el espectro de interferencia promedio @@ -1743,7 +1819,6 @@ class removeInterference(Operation): if (cnoiseid > 0): jspc_interf[noiseid] = 0 - # Expandiendo los perfiles a limpiar if (cinterfid > 0): new_interfid = ( @@ -1757,16 +1832,14 @@ class removeInterference(Operation): for ip in range(new_cinterfid): ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort() - jspc_interf[new_interfid[ip] - ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]] + jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]] - jspectra[ich, :, ind_hei] = jspectra[ich, :, - ind_hei] - jspc_interf # Corregir indices + jspectra[ich, :, ind_hei] = jspectra[ich, :,ind_hei] - jspc_interf # Corregir indices # Removiendo la interferencia del punto de mayor interferencia ListAux = jspc_interf[mask_prof].tolist() maxid = ListAux.index(max(ListAux)) - + print(cinterfid) if cinterfid > 0: for ip in range(cinterfid * (interf == 2) - 1): ind = (jspectra[ich, interfid[ip], :] < tmp_noise * @@ -1783,16 +1856,15 @@ class removeInterference(Operation): for id1 in range(4): xx[:, id1] = ind[id1]**numpy.asarray(list(range(4))) - xx_inv = numpy.linalg.inv(xx) xx = xx_inv[:, 0] ind = (ind + maxid + num_mask_prof) % num_mask_prof yy = jspectra[ich, mask_prof[ind], :] - jspectra[ich, mask_prof[maxid], :] = numpy.dot( - yy.transpose(), xx) + jspectra[ich, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx) indAux = (jspectra[ich, :, :] < tmp_noise * (1 - 1 / numpy.sqrt(num_incoh))).nonzero() + print(indAux) jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \ (1 - 1 / numpy.sqrt(num_incoh)) @@ -1856,7 +1928,7 @@ class removeInterference(Operation): return 1 - def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1): + def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1, minHei=None, maxHei=None): self.dataOut = dataOut @@ -1864,7 +1936,8 @@ class removeInterference(Operation): self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None) elif mode == 2: self.removeInterference2() - + elif mode == 3: + self.removeInterference3(min_hei=minHei, max_hei=maxHei) return self.dataOut @@ -2051,7 +2124,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 a22b44f..400fbec 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 numpy #import copy from schainpy.model.data import _noise @@ -413,6 +413,42 @@ class printAttribute(Operation): if hasattr(dataOut, attr): log.log(getattr(dataOut, attr), attr) +class cleanHeightsInterf(Operation): + + def run(self, dataOut, heightsList, repeats=0, step=0, factor=1, idate=None, startH=None, endH=None): + + #print(dataOut.data.shape) + + startTime = datetime.datetime.combine(idate,startH) + endTime = datetime.datetime.combine(idate,endH) + currentTime = datetime.datetime.fromtimestamp(dataOut.utctime) + + if currentTime < startTime or currentTime > endTime: + return dataOut + + wMask = numpy.asarray(factor) + wMask = numpy.tile(wMask,(repeats+2)) + #print(wMask) + heights = [float(hei) for hei in heightsList] + for r in range(repeats): + heights += [ (h+(step*(r+1))) for h in heights] + #print(heights) + heiList = dataOut.heightList + heights_indx = [getHei_index(h,h,heiList)[0] for h in heights] + + for ch in range(dataOut.data.shape[0]): + i = 0 + for hei in heights_indx: + #print(dataOut.data[ch,hei]) + if dataOut.data.ndim < 3: + dataOut.data[ch,hei-1] = (dataOut.data[ch,hei-1])*wMask[i] + else: + dataOut.data[ch,:,hei-1] = (dataOut.data[ch,:,hei-1])*wMask[i] + #print("done") + i += 1 + + return dataOut + class interpolateHeights(Operation): @@ -1785,7 +1821,7 @@ class SSheightProfiles2(Operation): residue = (self.__nHeis - self.nsamples) % self.step 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)) + print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,self.__nProfiles - self.nsamples,residue)) self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] self.init_range = dataOut.heightList[0] @@ -1808,14 +1844,15 @@ class SSheightProfiles2(Operation): if repeat is not None: code_block = numpy.repeat(code_block, repeats=repeat, axis=1) - 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 - else: - self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:] + if data.ndim < 3: + data = data.reshape(self.__nChannels,1,self.__nHeis ) + #print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape, code_block.shape) + for ch in range(self.__nChannels): + for i in range(int(self.new_nHeights)): #nuevas alturas + if code is not None: + self.buffer[ch,i,:] = data[ch,:,i*self.step:i*self.step + self.nsamples]*code_block + else: + self.buffer[ch,i,:] = data[ch,:,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,:,:]) @@ -1824,7 +1861,7 @@ class SSheightProfiles2(Operation): def run(self, dataOut, step, nsamples, code = None, repeat = None): - + # print("running") if dataOut.flagNoData == True: return dataOut dataOut.flagNoData = True @@ -1944,6 +1981,7 @@ class removeProfileByFaradayHS(Operation): self.nChannels = dataOut.nChannels self.nHeights = dataOut.nHeights + self.test_counter = 0 def filterSatsProfiles(self): data = self.__buffer_data @@ -2020,6 +2058,137 @@ class removeProfileByFaradayHS(Operation): self.outliers_IDs_list = numpy.unique(outlier_loc_index) return data + def cleanSpikesFFT2D(self): + incoh_int = 10 + norm_img = 75 + import matplotlib.pyplot as plt + import datetime + import cv2 + data = self.__buffer_data + print("cleaning shape inpt: ",data.shape) + self.__buffer_data = [] + + + channels , profiles, heights = data.shape + len_split_prof = profiles / incoh_int + + + for ch in range(channels): + data_10 = numpy.split(data[ch, :, self.minHei_idx:], incoh_int, axis=0) # división de los perfiles + print("splited data: ",len(data_10)," -> ", data_10[0].shape) + int_img = None + i_count = 0 + n_x, n_y = data_10[0].shape + for s_data in data_10: #porciones de espectro + spectrum = numpy.fft.fft2(s_data, axes=(0,1)) + z = numpy.abs(spectrum) + mg = z[2:n_y,:] #omitir dc y adjunto + dat = numpy.log10(mg.T) + i_count += 1 + if i_count == 1: + int_img = dat + else: + int_img += dat + #print(i_count) + + min, max = int_img.min(), int_img.max() + int_img = ((int_img-min)*255/(max-min)).astype(numpy.uint8) + + cv2.imshow('integrated image', int_img) #numpy.fft.fftshift(img)) + cv2.waitKey(0) + ##################################################################### + kernel_h = numpy.zeros((3,3)) # + kernel_h[0, :] = -2 + kernel_h[1, :] = 3 + kernel_h[2, :] = -2 + + + kernel_5h = numpy.zeros((5,5)) # + kernel_5h[0, :] = -2 + kernel_5h[1, :] = -1 + kernel_5h[2, :] = 5 + kernel_5h[3, :] = -1 + kernel_5h[4, :] = -2 + + ##################################################################### + sharp_img = cv2.filter2D(src=int_img, ddepth=-1, kernel=kernel_5h) + # cv2.imshow('sharp image h ', sharp_img) + # cv2.waitKey(0) + ##################################################################### + horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5,1)) #11 + ##################################################################### + detected_lines_h = cv2.morphologyEx(sharp_img, cv2.MORPH_OPEN, horizontal_kernel, iterations=1) + # cv2.imshow('lines horizontal', detected_lines_h) #numpy.fft.fftshift(detected_lines_h)) + # cv2.waitKey(0) + ##################################################################### + ret, detected_lines_h = cv2.threshold(detected_lines_h, 200, 255, cv2.THRESH_BINARY)# + cv2.imshow('binary img', detected_lines_h) #numpy.fft.fftshift(detected_lines_h)) + cv2.waitKey(0) + ##################################################################### + cnts_h, h0 = cv2.findContours(detected_lines_h, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) + ##################################################################### + h_line_index = [] + v_line_index = [] + + #cnts_h += cnts_h_s #combine large and small lines + + # line indexes x1, x2, y + for c in cnts_h: + #print(c) + if len(c) < 3: #contorno linea + x1 = c[0][0][0] + x2 = c[1][0][0] + if x1 > 5 and x2 < (n_x-5) : + start = incoh_int + (x1 * incoh_int) + end = incoh_int + (x2 * incoh_int) + h_line_index.append( [start, end, c[0][0][1]] ) + + #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1]) + else: #contorno poligono + pairs = numpy.asarray([c[n][0] for n in range(len(c))]) + y = numpy.unique(pairs[:,1]) + x = numpy.unique(pairs[:,0]) + #print(x) + for yk in y: + x0 = x[0] + if x0 < 8: + x0 = 10 + #print(x[0], x[-1], yk) + h_line_index.append( [x0, x[-1], yk]) + #print("x1, x2, y ->p ", x[0], x[-1], yk) + ################################################################### + #print("Cleaning") + # # clean Spectrum + spectrum = numpy.fft.fft2(data[ch,:,self.minHei_idx:], axes=(0,1)) + z = numpy.abs(spectrum) + phase = numpy.angle(spectrum) + print("Total Horizontal", len(h_line_index)) + if len(h_line_index) < 75 : + for x1, x2, y in h_line_index: + print(x1, x2, y) + z[x1:x2,y] = 0 + + + spcCleaned = z * numpy.exp(1j*phase) + + dat2 = numpy.log10(z[1:-1,:].T) + min, max =dat2.min(), dat2.max() + print(min, max) + img2 = ((dat2-min)*255/(max-min)).astype(numpy.uint8) + cv2.imshow('cleaned', img2) #numpy.fft.fftshift(img_cleaned)) + cv2.waitKey(0) + cv2.destroyAllWindows() + + data[ch,:,self.minHei_idx:] = numpy.fft.ifft2(spcCleaned, axes=(0,1)) + + + #print("cleanOutliersByBlock Done", data.shape) + self.__buffer_data = data + return data + + + + def cleanOutliersByBlock(self): import matplotlib.pyplot as plt import datetime @@ -2027,10 +2196,10 @@ class removeProfileByFaradayHS(Operation): #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[:,:,self.minHei_idx:], axes=(0,2)) + + spectrum = numpy.fft.fft2(data[:,:,self.minHei_idx:], axes=(1,2)) print("spc : ",spectrum.shape) (nch,nsamples, nh) = spectrum.shape data2 = None @@ -2045,134 +2214,340 @@ class removeProfileByFaradayHS(Operation): freqv = numpy.fft.fftfreq(nh, d=dt1) freqh = numpy.fft.fftfreq(self.n, d=dt2) - # freqv = numpy.fft.fftshift(numpy.fft.fftfreq(nh, d=dt1)) - # freqh = numpy.fft.fftshift(numpy.fft.fftfreq(self.n, d=dt2)) - # + z = numpy.abs(spectrum[ch,:,:]) + phase = numpy.angle(spectrum[ch,:,:]) + z1 = z[0,:] + #print("shape z: ", z.shape, nsamples) + dat = numpy.log10(z[1:nsamples,:].T) + pdat = numpy.log10(phase.T) + #print("dat mean",dat.mean()) + mean, min, max = dat.mean(), dat.min(), dat.max() + img = ((dat-min)*200/(max-min)).astype(numpy.uint8) + # print(img.shape) + cv2.imshow('image', img) #numpy.fft.fftshift(img)) + cv2.waitKey(0) - z = numpy.abs(spectrum[ch,:,:]) - phase = numpy.angle(spectrum[ch,:,:]) + ''' #FUNCIONA LINEAS PEQUEÑAS + kernel_5h = numpy.zeros((5,3)) # + kernel_5h[0, :] = 2 + kernel_5h[1, :] = 1 + kernel_5h[2, :] = 0 + kernel_5h[3, :] = -1 + kernel_5h[4, :] = -2 - print("shape z: ", z.shape) + sharp_imgh = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_5h) + cv2.imshow('sharp image h',sharp_imgh) + cv2.waitKey(0) + horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (20,1)) + detected_lines_h = cv2.morphologyEx(sharp_imgh, cv2.MORPH_OPEN, horizontal_kernel, iterations=1) + #detected_lines_h = cv2.medianBlur(detected_lines_h, 3) + #detected_lines_h = cv2.filter2D(src=img, ddepth=-1, kernel=kernel) + cv2.imshow('lines h gray', detected_lines_h) + cv2.waitKey(0) + reth, detected_lines_h = cv2.threshold(detected_lines_h, 90, 255, cv2.THRESH_BINARY) + cv2.imshow('lines h ', detected_lines_h) + cv2.waitKey(0) + ''' - # Find all peaks higher than the 98th percentile - peaks = z < numpy.percentile(z, 99) - #print(peaks) - # Set those peak coefficients to zero - #spectrum2 = spectrum * peaks.astype(int) - #data2 = numpy.fft.ifft2(spectrum2,axes=(0,2)) - dat = numpy.log10(z.T) - pdat = numpy.log10(phase.T) + ''' + kernel_3h = numpy.zeros((3,10)) #10 + kernel_3h[0, :] = -1 + kernel_3h[1, :] = 2 + kernel_3h[2, :] = -1 - min, max = dat.min(), dat.max() - #img = ((dat-min)*255/(max-min)).astype(numpy.uint8) - img = ((dat-min)*200/(max-min)).astype(numpy.uint8) - #img = 255 - ((dat-min)*255/max).astype(numpy.uint8) + kernel_h = numpy.zeros((3,20)) #20 + kernel_h[0, :] = -1 + kernel_h[1, :] = 2 + kernel_h[2, :] = -1 - #print(img.shape) - # cv2.imshow('image', numpy.fft.fftshift(img)) - # cv2.waitKey(0) + kernel_v = numpy.zeros((30,3)) #30 + kernel_v[:, 0] = -1 + kernel_v[:, 1] = 2 + kernel_v[:, 2] = -1 + + kernel_4h = numpy.zeros((4,20)) # + kernel_4h[0, :] = 1 + kernel_4h[1, :] = 0 + kernel_4h[2, :] = 0 + kernel_4h[3, :] = -1 - spikes = numpy.where(img > 230, True, False) - #print("spikes: ", spikes.sum()) + kernel_5h = numpy.zeros((5,30)) # + kernel_5h[0, :] = 2 + kernel_5h[1, :] = 1 + kernel_5h[2, :] = 0 + kernel_5h[3, :] = -1 + kernel_5h[4, :] = -2 - #img = cv2.medianBlur(img,3) - #cv2.imshow('image', cv2.resize(img.astype('float32'), (600, 600))) - kernel3 = numpy.zeros((3,15)) - kernel3[0, :] = -1 - kernel3[1, :] = 2 - kernel3[2, :] = -1 - sharp_img = cv2.filter2D(src=img, ddepth=-1, kernel=kernel3) + sharp_img0 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_3h) + # cv2.imshow('sharp image small h',sharp_img0) # numpy.fft.fftshift(sharp_img1)) + # cv2.waitKey(0) + + sharp_img1 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_h) + # cv2.imshow('sharp image h',sharp_img1) # numpy.fft.fftshift(sharp_img1)) + # cv2.waitKey(0) + sharp_img2 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_v) + # cv2.imshow('sharp image v', sharp_img2) #numpy.fft.fftshift(sharp_img2)) + # cv2.waitKey(0) - # cv2.imshow('sharp image', numpy.fft.fftshift(sharp_img)) + sharp_imgw = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_4h) + # cv2.imshow('sharp image h wide', sharp_imgw) #numpy.fft.fftshift(sharp_img2)) # cv2.waitKey(0) - ret, thresh = cv2.threshold(sharp_img, 127, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU) - # #thresh = cv2.adaptiveThreshold(sharp_img,255,cv2.ADAPTIVE_THRESH_MEAN_C,cv2.THRESH_BINARY_INV,21,0) - # cv2.imshow('binary img', numpy.fft.fftshift(thresh)) + sharp_imgwl = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_5h, borderType = cv2.BORDER_ISOLATED) + cv2.imshow('sharp image h long wide', sharp_imgwl) #numpy.fft.fftshift(sharp_img2)) + cv2.waitKey(0) + + # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/sharp_h.jpg', sharp_img1) + # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/sharp_v.jpg', sharp_img2) + # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/input_img.jpg', img) + + ########################small horizontal + horizontal_kernel_s = cv2.getStructuringElement(cv2.MORPH_RECT, (11,1)) #11 + ######################## horizontal + horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (30,1)) #30 + ######################## vertical + vertical_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (1,50)) #50 + ######################## horizontal wide + horizontal_kernel_w = cv2.getStructuringElement(cv2.MORPH_RECT, (30,1)) # 30 + + horizontal_kernel_expand = cv2.getStructuringElement(cv2.MORPH_RECT, (3,3)) # + + horizontal_kernel_wl = cv2.getStructuringElement(cv2.MORPH_RECT, (50,1)) # + + detected_lines_h_s = cv2.morphologyEx(sharp_img0, cv2.MORPH_OPEN, horizontal_kernel_s, iterations=7) #7 + detected_lines_h = cv2.morphologyEx(sharp_img1, cv2.MORPH_OPEN, horizontal_kernel, iterations=7) #7 + detected_lines_v = cv2.morphologyEx(sharp_img2, cv2.MORPH_OPEN, vertical_kernel, iterations=7) #7 + detected_lines_h_w = cv2.morphologyEx(sharp_imgw, cv2.MORPH_OPEN, horizontal_kernel_w, iterations=5) #5 + + detected_lines_h_wl = cv2.morphologyEx(sharp_imgwl, cv2.MORPH_OPEN, horizontal_kernel_wl, iterations=5) # + detected_lines_h_wl = cv2.filter2D(src=detected_lines_h_wl, ddepth=-1, kernel=horizontal_kernel_expand) + + # cv2.imshow('lines h small gray', detected_lines_h_s) #numpy.fft.fftshift(detected_lines_h)) + # cv2.waitKey(0) + # cv2.imshow('lines h gray', detected_lines_h) #numpy.fft.fftshift(detected_lines_h)) # cv2.waitKey(0) - # #Remove horizontal - horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (10,1)) - detected_lines = cv2.morphologyEx(sharp_img, cv2.MORPH_OPEN, horizontal_kernel, iterations=8) - cnts, h = cv2.findContours(detected_lines, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) + # cv2.imshow('lines v gray', detected_lines_v) #numpy.fft.fftshift(detected_lines_h)) + # cv2.waitKey(0) + # cv2.imshow('lines h wide gray', detected_lines_h_w) #numpy.fft.fftshift(detected_lines_h)) + # cv2.waitKey(0) + cv2.imshow('lines h long wide gray', detected_lines_h_wl) #numpy.fft.fftshift(detected_lines_h)) + cv2.waitKey(0) + + reth_s, detected_lines_h_s = cv2.threshold(detected_lines_h_s, 85, 255, cv2.THRESH_BINARY)# 85 + reth, detected_lines_h = cv2.threshold(detected_lines_h, 30, 255, cv2.THRESH_BINARY) #30 + retv, detected_lines_v = cv2.threshold(detected_lines_v, 30, 255, cv2.THRESH_BINARY) #30 + reth_w, detected_lines_h_w = cv2.threshold(detected_lines_h_w, 35, 255, cv2.THRESH_BINARY)# + reth_wl, detected_lines_h_wl = cv2.threshold(detected_lines_h_wl, 200, 255, cv2.THRESH_BINARY)# + + cnts_h_s, h0 = cv2.findContours(detected_lines_h_s, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) + cnts_h, h1 = cv2.findContours(detected_lines_h, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) + cnts_v, h2 = cv2.findContours(detected_lines_v, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) + cnts_h_w, h3 = cv2.findContours(detected_lines_h_w, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) + cnts_h_wl, h4 = cv2.findContours(detected_lines_h_wl, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE) + #print("horizontal ", cnts_h) + #print("vertical ", cnts_v) + # cnts, h = cv2.findContours(detected_lines_h, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE) # print(cnts) - # cv2.imshow('lines ', numpy.fft.fftshift(detected_lines)) + # cv2.imshow('lines h wide', detected_lines_h_w) #numpy.fft.fftshift(detected_lines_h)) + # cv2.waitKey(0) + cv2.imshow('lines h wide long ', detected_lines_h_wl) #numpy.fft.fftshift(detected_lines_v)) + cv2.waitKey(0) + # cv2.imshow('lines h small', detected_lines_h_s) #numpy.fft.fftshift(detected_lines_h)) + # cv2.waitKey(0) + # cv2.imshow('lines h ', detected_lines_h) #numpy.fft.fftshift(detected_lines_h)) # cv2.waitKey(0) + # cv2.imshow('lines v ', detected_lines_v) #numpy.fft.fftshift(detected_lines_v)) + # cv2.waitKey(0) + + # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/lines_h.jpg', detected_lines_h) + # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/lines_v.jpg', detected_lines_v) + #cnts = cnts[0] if len(cnts) == 2 else cnts[1] - y_line_index = numpy.asarray([ [c[0][0][0],c[1][0][0], c[0][0][1]] for c in cnts ]) - - # print(y_line_index) - u, index, counts = numpy.unique(y_line_index[:,2], return_index=True, return_counts=True) - print(u, index, counts) - readed_lines = [] - for x1, x2, y, YY, n in zip(y_line_index, u, counts): - if y in readed_lines: - pass - if n > 1: - len_line = self.n - else: - len_line = len(z[x1:x2,y]) - #print(y, nh, self.n) - if y != (nh-1): - list = [ ((z[n, y-1] + z[n,y+1])/2) for n in range(len_line)] - else: - list = [ ((z[n, y-1] + z[n,0])/2) for n in range(len_line)] - if n > 1: - z[:,y] = numpy.asarray(list) - else: - z[x1:x2,y] = numpy.asarray(list) + #y_line_index = numpy.asarray([ [c[0][0][0],c[1][0][0], c[0][0][1]] for c in cnts_v ]) + h_line_index = [] + v_line_index = [] - readed_lines.append(y) + cnts_h += cnts_h_s #combine large and small lines + + # line indexes x1, x2, y + for c in cnts_h: + #print(c) + if len(c) < 3: #contorno linea + x1 = c[0][0][0] + if x1 < 8: + x1 = 10 + h_line_index.append( [x1,c[1][0][0], c[0][0][1]] ) + #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1]) + else: #contorno poligono + pairs = numpy.asarray([c[n][0] for n in range(len(c))]) + y = numpy.unique(pairs[:,1]) + x = numpy.unique(pairs[:,0]) + #print(x) + for yk in y: + x0 = x[0] + if x0 < 8: + x0 = 10 + #print(x[0], x[-1], yk) + h_line_index.append( [x0, x[-1], yk]) + #print("x1, x2, y ->p ", x[0], x[-1], yk) + for c in cnts_h_w: + #print(c) + if len(c) < 3: #contorno linea + x1 = c[0][0][0] + if x1 < 8: + x1 = 10 + y = c[0][0][1] - 2 # se incrementa 2 líneas x el filtro + h_line_index.append( [x1,c[1][0][0],y] ) + #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1]) + else: #contorno poligono + pairs = numpy.asarray([c[n][0] for n in range(len(c))]) + y = numpy.unique(pairs[:,1]) + x = numpy.unique(pairs[:,0]) + #print(x) + for yk in y: + + x0 = x[0] + if x0 < 8: + x0 = 10 + h_line_index.append( [x0, x[-1], yk-2]) + + for c in cnts_h_wl: # # revisar + #print(c) + if len(c) < 3: #contorno linea + x1 = c[0][0][0] + if x1 < 8: + x1 = 10 + y = c[0][0][1] - 2 # se incrementa 2 líneas x el filtro + h_line_index.append( [x1,c[1][0][0],y] ) + #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1]) + else: #contorno poligono + pairs = numpy.asarray([c[n][0] for n in range(len(c))]) + y = numpy.unique(pairs[:,1]) + x = numpy.unique(pairs[:,0]) + for yk in range(y[-1]-y[0]): + y_k = yk +y[0] + + x0 = x[0] + if x0 < 8: + x0 = 10 + h_line_index.append( [x0, x[-1], y_k-2]) + + print([[c[0][0][1],c[1][0][1], c[0][0][0] ] for c in cnts_v]) + # line indexes y1, y2, x + for c in cnts_v: + if len(c) < 3: #contorno linea + v_line_index.append( [c[0][0][1],c[1][0][1], c[0][0][0] ] ) + else: #contorno poligono + pairs = numpy.asarray([c[n][0] for n in range(len(c))]) + #print(pairs) + y = numpy.unique(pairs[:,1]) + x = numpy.unique(pairs[:,0]) + + for xk in x: + #print(x[0], x[-1], yk) + v_line_index.append( [y[0],y[-1], xk]) + + ################################################################### + # # clean Horizontal + print("Total Horizontal", len(h_line_index)) + if len(h_line_index) < 75 : + for x1, x2, y in h_line_index: + #print("cleaning ",x1, x2, y) + len_line = x2 - x1 + if y > 10 and y < (nh -10): + # if y != (nh-1): + # list = [ ((z[n, y-1] + z[n,y+1])/2) for n in range(len_line)] + # else: + # list = [ ((z[n, y-1] + z[n,0])/2) for n in range(len_line)] + # + # z[x1:x2,y] = numpy.asarray(list) + z[x1-5:x2+5,y:y+1] = 0 + + # clean vertical + for y1, y2, x in v_line_index: + len_line = y2 - y1 + #print(x) + if x > 0 and x < (nsamples-2): + # if x != (nsamples-1): + # list = [ ((z[x-2, n] + z[x+2,n])/2) for n in range(len_line)] + # else: + # list = [ ((z[x-2, n] + z[1,n])/2) for n in range(len_line)] + # + # #z[x-1:x+1,y1:y2] = numpy.asarray(list) + # + z[x+1,y1:y2] = 0 + + ''' + #z[: ,[215, 217, 221, 223, 225, 340, 342, 346, 348, 350, 465, 467, 471, 473, 475]]=0 + z[1: ,[112, 114, 118, 120, 122, 237, 239, 245, 247, 249, 362, 364, 368, 370, 372]]=0 + # z[: ,217]=0 + # z[: ,221]=0 + # z[: ,223]=0 + # z[: ,225]=0 dat2 = numpy.log10(z.T) - img = ((dat2-min)*255/(max-min)).astype(numpy.uint8) + #print(dat2) + max = dat2.max() + #print(" min, max ", max, min) + img2 = ((dat2-min)*255/(max-min)).astype(numpy.uint8) + #img_cleaned = img2.copy() + #cv2.drawContours(img2, cnts_h, -1, (255,255,255), 1) + #cv2.drawContours(img2, cnts_v, -1, (255,255,255), 1) - for c in cnts: - #print(c) - cv2.drawContours(img, [c], -1, (255,255,255), 2) - #print("C: ",c[0][1]) spcCleaned = z * numpy.exp(1j*phase) #print(spcCleaned) - cv2.imshow('image2', numpy.fft.fftshift(img)) + # cv2.imshow('image contours', img2) #numpy.fft.fftshift(img)) + # cv2.waitKey(0) + + cv2.imshow('cleaned', img2) #numpy.fft.fftshift(img_cleaned)) cv2.waitKey(0) + # # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/cleaned_{}.jpg'.format(self.test_counter), img2) cv2.destroyAllWindows() + # self.test_counter += 1 - m = numpy.mean(dat) - o = numpy.std(dat) - fig, ax = plt.subplots(1,2,figsize=(12, 6)) - #X, Y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv)) - X, Y = numpy.meshgrid(numpy.fft.fftshift(freqh),numpy.fft.fftshift(freqv)) + #print("DC difference " ,z1 - z[0,:]) - colormap = 'jet' - #c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o)) - c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = 4.2, vmax = 5.0) - fig.colorbar(c, ax=ax[0]) - - - #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) - #print("aqui estoy2",dat2[:,:,0].shape) - c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(dat2), cmap =colormap, vmin = 4.2, vmax = 5.0) - #c = ax[1].pcolormesh(x, y, dat2[:,:,0], cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o)-1) - #print("aqui estoy3") - fig.colorbar(c, ax=ax[1]) - plt.show() + # m = numpy.mean(dat) + # o = numpy.std(dat) + # print("mean ", m, " std ", o) + # fig, ax = plt.subplots(1,2,figsize=(12, 6)) + # #X, Y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv)) + # X, Y = numpy.meshgrid(numpy.fft.fftshift(freqh),numpy.fft.fftshift(freqv)) + # + # colormap = 'jet' + # #c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o)) + # #c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = 6.5, vmax = 6.8) + # c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = (m-2*o), vmax = (m+1.5*o)) + # fig.colorbar(c, ax=ax[0]) + # + # + # #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) + # #print("aqui estoy2",dat2[:,:,0].shape) + # #c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(pdat), cmap =colormap, vmin = 4.2, vmax = 5.0) + # c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat2), cmap =colormap, vmin = (m-2*o), vmax = (m+1.5*o)) + # #c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(pdat), cmap =colormap ) #, vmin = 0.0, vmax = 0.5) + # #c = ax[1].pcolormesh(x, y, dat2[:,:,0], cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o)-1) + # #print("aqui estoy3") + # fig.colorbar(c, ax=ax[1]) + # plt.show() spectrum[ch,:,:] = spcCleaned @@ -2180,10 +2555,10 @@ class removeProfileByFaradayHS(Operation): - data[:,:,self.minHei_idx:] = numpy.fft.ifft2(spectrum, axes=(0,2)) + data[:,:,self.minHei_idx:] = numpy.fft.ifft2(spectrum, axes=(1,2)) #print("cleanOutliersByBlock Done", data.shape) - + self.__buffer_data = data return data @@ -2212,6 +2587,7 @@ class removeProfileByFaradayHS(Operation): if self.__profIndex == self.n: #print("apnd : ",data) #dataBlock = self.cleanOutliersByBlock() + #dataBlock = self.cleanSpikesFFT2D() dataBlock = self.filterSatsProfiles() self.__dataReady = True @@ -2288,7 +2664,7 @@ class removeProfileByFaradayHS(Operation): 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, dataOut.flagNoData) self.n_prof_released += 1 if self.end_prof >= (self.n +self.lenProfileOut):