schain Commit - r1655:13a9a26f11d7 · Jicamarca Repository

upgrade filterbyHeights, fix initial parameters. removeProfileSats2 improved debug

joabAM -

r1655:13a9a26f11d7

parent child

Context file:

r1655:13a9a26f11d7

Collapse all files

schainpy/model/graphics/jroplot_base.py +4 -2

                     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.clear_figures()
-                        self.__plot()
                     self.__update(dataOut, tm)

schainpy/model/proc/jroproc_voltage.py +67 -41

             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
+                        window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / self.deltaHeight
-                    newdelta = deltaHeight * window
+                    if not self.ifConfig: #and dataOut.useInputBuffer:
-                    r = dataOut.nHeights % window
+                        self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
-                    newheights = (dataOut.nHeights-r)/window
+                        self.ifConfig = True
+                        self.newdelta = self.deltaHeight * window
-                    if newheights <= 1:
+                        self.r = dataOut.nHeights % window
-                        raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(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 = dataOut.data[:, :, 0:int(self.nHeights-self.r)]
-                        buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
+                        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 = dataOut.data[:,0:int(self.nHeights-self.r)]
-                        buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
+                        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):
                 __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):
                     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 ):
                     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!")
                         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)
                         #print(self.min_ref,self.max_ref)
                         import scipy.signal
-                        b, a = scipy.signal.butter(3, 0.1)
+                        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])).real
+                        #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
-                        power = (data[c,:,self.minHei_idx:self.maxHei_idx] * numpy.conjugate(data[c,:,self.minHei_idx:self.maxHei_idx])).real
+                        #power = (data[c,:,self.minHei_idx:self.maxHei_idx] * numpy.conjugate(data[c,:,self.minHei_idx:self.maxHei_idx]))
-                        heis = len(power[0,:])
+                        power = numpy.abs(data[c,:,self.minHei_idx:self.maxHei_idx])
-                        power = power.mean(axis=1)
+                        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())
                             # plt.grid()
                             # plt.show()
-                        # fig, ax = plt.subplots()
+                        if self.debug:
-                        # ax.plot(fpower)
+                            fig, ax = plt.subplots()
-                        # ax.axhline(th, color='r')
+                            ax.plot(fpower, label="power")
-                        # ax.axhline(std, color='b')
+                            #ax.plot(fnoise, label="noise ref")
-                        # plt.grid()
+                            ax.axhline(th, color='g', label="th")
-                        # plt.show()
+                            #ax.axhline(std, color='b', label="mean")
+                            ax.legend()
+                            plt.grid()
+                            plt.show()
                         #print(indexes)
                         indexesYagi = indexesYagi[ (indexesYagi >= 0) & (indexesYagi<profiles)]
                         indexesYagi = numpy.unique(indexesYagi)
+                    #print("indexes: " ,indexes)
                     outs_lines = numpy.unique(indexes)
+                    #print(outs_lines)
                     #Agrupando el histograma de outliers,
                     my_bins = numpy.linspace(0,int(profiles), div, endpoint=True)
                     hist, bins =  numpy.histogram(outs_lines,bins=my_bins)
-                    hist_outliers_indexes = numpy.where(hist >= self.thHistOutlier)       #es outlier
+                    #print("hist: ",hist)
-                    hist_outliers_indexes = hist_outliers_indexes[0]
+                    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)
                                 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')

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages