From 13a9a26f11d7ad38cf35648b1532e522badb4782 2023-07-12 16:23:03 From: joabAM Date: 2023-07-12 16:23:03 Subject: [PATCH] upgrade filterbyHeights, fix initial parameters. removeProfileSats2 improved debug --- diff --git a/schainpy/model/graphics/jroplot_base.py b/schainpy/model/graphics/jroplot_base.py index 2860f77..d157fb2 100644 --- a/schainpy/model/graphics/jroplot_base.py +++ b/schainpy/model/graphics/jroplot_base.py @@ -665,11 +665,13 @@ class Plot(Operation): tm = getattr(dataOut, self.attr_time) if self.data and 'time' in self.xaxis and (tm - self.tmin) >= self.xrange*60*60: + self.clear_figures() self.save_time = tm + self.__plot() self.tmin += self.xrange*60*60 self.data.setup() - self.clear_figures() - self.__plot() + #self.clear_figures() + self.__update(dataOut, tm) diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 33979d5..8bf0b62 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -320,45 +320,57 @@ class selectHeights(Operation): class filterByHeights(Operation): - + ifConfig=False + deltaHeight = None + newdelta=None + newheights=None + r=None + h0=None + nHeights=None def run(self, dataOut, window): - - deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] - + + + # print("1",dataOut.data.shape) + # print(dataOut.nHeights) if window == None: - window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight - - newdelta = deltaHeight * window - r = dataOut.nHeights % window - newheights = (dataOut.nHeights-r)/window - - if newheights <= 1: - raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window)) + window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / self.deltaHeight + + if not self.ifConfig: #and dataOut.useInputBuffer: + self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] + self.ifConfig = True + self.newdelta = self.deltaHeight * window + self.r = dataOut.nHeights % window + self.newheights = (dataOut.nHeights-self.r)/window + self.h0 = dataOut.heightList[0] + self.nHeights = dataOut.nHeights + if self.newheights <= 1: + raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window)) if dataOut.flagDataAsBlock: """ Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis] """ - buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)] - buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window) + buffer = dataOut.data[:, :, 0:int(self.nHeights-self.r)] + buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(self.nHeights/window), window) buffer = numpy.sum(buffer,3) else: - buffer = dataOut.data[:,0:int(dataOut.nHeights-r)] - buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window)) + buffer = dataOut.data[:,0:int(self.nHeights-self.r)] + buffer = buffer.reshape(dataOut.nChannels,int(self.nHeights/window),int(window)) buffer = numpy.sum(buffer,2) dataOut.data = buffer - dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta + dataOut.heightList = self.h0 + numpy.arange( self.newheights )*self.newdelta dataOut.windowOfFilter = window #update Processing Header: dataOut.processingHeaderObj.heightList = dataOut.heightList dataOut.processingHeaderObj.nWindows = window - + return dataOut + class setH0(Operation): def run(self, dataOut, h0, deltaHeight = None): @@ -3276,7 +3288,7 @@ class RemoveProfileSats2(Operation): __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights', - 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels', + 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels','cohFactor', '__count_exec','__initime','__dataReady','__ipp', 'minRef', 'maxRef', 'debug','prev_pnoise') def __init__(self, **kwargs): @@ -3284,7 +3296,7 @@ class RemoveProfileSats2(Operation): self.isConfig = False self.currentTime = None - def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=15,minHei=None, maxHei=None, nBins=10, + def setup(self,dataOut, n=None , navg=0.9, profileMargin=50,thHistOutlier=15,minHei=None, maxHei=None, nBins=10, minRef=None, maxRef=None, debug=False, remYagi=False, nProfYagi = 0, offYagi=0, minHJULIA=None, maxHJULIA=None, idate=None,startH=None,endH=None ): @@ -3323,7 +3335,7 @@ class RemoveProfileSats2(Operation): self.test_counter = 0 self.debug = debug self.remYagi = remYagi - + self.cohFactor = dataOut.nCohInt if self.remYagi : if minHJULIA==None or maxHJULIA==None: raise ValueError("Parameters minHYagi and minHYagi are necessary!") @@ -3338,7 +3350,7 @@ class RemoveProfileSats2(Operation): self.startTime = datetime.datetime.combine(idate,startH) self.endTime = datetime.datetime.combine(idate,endH) - log.warning("Be careful with the selection of parameters for sat removal! Is recommended to \ + log.warning("Be careful with the selection of parameters for sats removal! It is avisable to \ activate the debug parameter in this operation for calibration", self.name) @@ -3370,14 +3382,19 @@ activate the debug parameter in this operation for calibration", self.name) #print(self.min_ref,self.max_ref) import scipy.signal - b, a = scipy.signal.butter(3, 0.1) - noise_ref = (data[c,:,self.min_ref:self.max_ref] * numpy.conjugate(data[c,:,self.min_ref:self.max_ref])).real + b, a = scipy.signal.butter(3, 0.5) + #noise_ref = (data[c,:,self.min_ref:self.max_ref] * numpy.conjugate(data[c,:,self.min_ref:self.max_ref])) + noise_ref = numpy.abs(data[c,:,self.min_ref:self.max_ref]) + lnoise = len(noise_ref[0,:]) + #print(noise_ref.shape) noise_ref = noise_ref.mean(axis=1) + #fnoise = noise_ref fnoise = scipy.signal.filtfilt(b, a, noise_ref) #noise_refdB = 10* numpy.log10(noise_ref) #print("Noise ",numpy.percentile(noise_ref,95)) - p85 = numpy.percentile(fnoise,83) + p95 = numpy.percentile(fnoise,90) mean_noise = fnoise.mean() + if self.prev_pnoise != None: if mean_noise < (1.5 * self.prev_pnoise) : self.prev_pnoise = mean_noise @@ -3390,20 +3407,24 @@ activate the debug parameter in this operation for calibration", self.name) - power = (data[c,:,self.minHei_idx:self.maxHei_idx] * numpy.conjugate(data[c,:,self.minHei_idx:self.maxHei_idx])).real - heis = len(power[0,:]) - power = power.mean(axis=1) + #power = (data[c,:,self.minHei_idx:self.maxHei_idx] * numpy.conjugate(data[c,:,self.minHei_idx:self.maxHei_idx])) + power = numpy.abs(data[c,:,self.minHei_idx:self.maxHei_idx]) + npower = len(power[0,:]) + #print(power.shape) + power = power.mean(axis=1) fpower = scipy.signal.filtfilt(b, a, power) #print(power.shape) #powerdB = 10* numpy.log10(power) - th = p85 + th = p95 #* 1.1 + #th = mean_noise index = numpy.where(fpower > th ) - #print("Noise ",mean_noise, p85) + #print("Noise ",mean_noise, p95) #print(index) - if index[0].size < int(self.navg*profiles): + + if index[0].size <= int(self.navg*profiles): indexes = numpy.append(indexes, index[0]) #print("sdas ", noise_ref.mean()) @@ -3423,12 +3444,15 @@ activate the debug parameter in this operation for calibration", self.name) # plt.grid() # plt.show() - # fig, ax = plt.subplots() - # ax.plot(fpower) - # ax.axhline(th, color='r') - # ax.axhline(std, color='b') - # plt.grid() - # plt.show() + if self.debug: + fig, ax = plt.subplots() + ax.plot(fpower, label="power") + #ax.plot(fnoise, label="noise ref") + ax.axhline(th, color='g', label="th") + #ax.axhline(std, color='b', label="mean") + ax.legend() + plt.grid() + plt.show() #print(indexes) @@ -3464,16 +3488,18 @@ activate the debug parameter in this operation for calibration", self.name) indexesYagi = indexesYagi[ (indexesYagi >= 0) & (indexesYagi= self.thHistOutlier) #es outlier - hist_outliers_indexes = hist_outliers_indexes[0] + #print("hist: ",hist) + hist_outliers_indexes = numpy.where(hist >= self.thHistOutlier)[0] #es outlier + #print(hist_outliers_indexes) if len(hist_outliers_indexes>0): hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1) @@ -3516,7 +3542,7 @@ activate the debug parameter in this operation for calibration", self.name) if self.remYagi : ax[0].vlines(indexesYagi,750,850, linestyles='dashed', label = 'yagi', color='m') else: - c = ax[0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70) + c = ax[0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = (70+2*self.cohFactor)) ax[0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w') #fig.colorbar(c) ax[0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')