Changeset 1161 for SHX


Ignore:
Timestamp:
05.02.2016 21:35:37 (4 years ago)
Author:
klaus
Message:

corrpick now working

Location:
SHX/trunk/SeismicHandler/commands
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • SHX/trunk/SeismicHandler/commands/corrpick.py

    r1160 r1161  
    66 
    77""" 
    8 Correlation picker. 
     8Smart correlation picker. Computes correlation of the selected time window with 
     9all other traces. Takes the N (N=3) largest correlation values and takes the 
     10best fitting a plane wave. If residuals stay over a tolerance value, no pick 
     11is set on this trace. 
    912""" 
    1013 
     
    3134    """ 
    3235    # Define number and type of parameters, possible qualifiers. 
    33     numberOfParameters = [4] 
     36    numberOfParameters = [4,5] 
    3437    parameterQueries = [ 
    3538        { 
     
    5558    ] 
    5659    known_qualifiers = [ 
     60        "CORRMAXIMA", 
     61        "MAXTOL", 
     62        "MAXLOOP", 
    5763    ] 
    5864    redraw = True 
     
    6773    #@timeit 
    6874    def run(self): 
    69         "Correlation picker." 
     75        "Smart correlation picker." 
    7076 
    7177        nummax = 3 
     78        if self.qualifiers["CORRMAXIMA"]: 
     79            nummax = int( self.qualifers["CORRMAXIMA"] ) 
     80        maxtol = 2. 
     81        if self.qualifiers["MAXTOL"]: 
     82            maxtol = float( self.qualifers["MAXTOL"] ) 
     83        maxloop = 100 
     84        if self.qualifiers["MAXLOOP"]: 
     85            maxloop = int( self.qualifers["MAXLOOP"] ) 
    7286        self.modifdict = {} 
    7387        self.maxmodif = 2*nummax 
     
    7892        self.phaselist = PhaseList() 
    7993        trace = traces[0] 
     94        # compute time windows for correlation 
    8095        _start = float( self.parameters[1] ) 
    8196        _end = float( self.parameters[2] ) 
    8297        _shift = float( self.parameters[3] ) 
     98        if len(self.parameters) > 4: 
     99            outvar = self.parameters[4][1:] 
     100        else: 
     101            outvar = None 
    83102        wdwlen = _end - _start 
    84103        trace_sname = "%s.%s" % (trace.stats.network,trace.stats.station) 
    85         # find phase 
     104        # find phase in selected window 
    86105        absstart = trace.stats.starttime + _start - trace.get_info('t-origin') 
    87106        absend = absstart + wdwlen 
     
    96115            shiftend = trace.stats.endtime 
    97116        _shift = ((shiftend - shiftstart) - wdwlen)/2. 
     117        # Prepare traces ad correlate. 
     118        # Results are in corrinfo, which is a dictionary on station names. 
     119        # Each station has a list of pairs (reltime,relcorrvalue), reltime is 
     120        # the time of the correlation maxima, relcorrvalue is the relative 
     121        # correlation coefficient (ratio to the largest coefficient). 
    98122        corrtrc = trace.slice( shiftstart, shiftend ) 
    99123        shiftlength = int( _shift / trace.stats.delta )/2 - 1 
     
    127151                except: 
    128152                    pass 
    129         corrinfo = self.findOptShifts( corrinfo ) 
     153        # Modify the maxima to fit a plane wave. 
     154        corrinfo, res = self.findOptShifts( corrinfo, maxtol, maxloop ) 
     155        # Define phases as given in the first list entry in each station. 
    130156        for sname in corrinfo.keys(): 
    131             reltime = corrinfo[sname][0][0] 
     157            if len(corrinfo[sname]) > 0: 
     158                reltime = corrinfo[sname][0][0] 
     159            else: 
     160                reltime = None 
    132161            self.phaselist.clearPhasenameStation( phase, sname, auto ) 
    133162            if reltime == None: 
     
    135164            self.phaselist.addPhasePick( sname, phasetime+reltime, name=phase, 
    136165                picktype=auto ) 
     166        if outvar: 
     167            self.symbols.set( outvar, "%g" % res ) 
    137168     
    138169    def findPhase( self, trc, start, end ): 
     170        "Determin name and time of phase in the selected time window." 
    139171        sname = "%s.%s" % (trc.stats.network,trc.stats.station) 
    140172        manual = self.phaselist.translatePicktype( 'manual' ) 
     
    150182        return (None,None) 
    151183     
    152     def findOptShifts( self, corrinfo ): 
     184    def findOptShifts( self, corrinfo, maxtol=2., maxloop=100 ): 
     185        """ 
     186        Modify the correlation maxima (cyclic exchange or delete) to match 
     187        a plane wave. 
     188        """ 
    153189        # find station locations 
    154190        optcorrinfo = None 
     
    168204            lats.append( lat ) 
    169205            lons.append( lon ) 
     206        # recompute to relative coos in km 
    170207        relx, rely = absToRelPositions( lons, lats ) 
    171         for i in range(50): 
     208        # Several tries to modify corrinfo. 
     209        nonecnt = 0 
     210        for i in range(maxloop): 
     211            # Collect rel. times from correlation maxima in the current order. 
     212            # 'None' values at stations with too large persistent residuals. 
    172213            reltimes = [] 
    173214            for s in names: 
     
    176217                else: 
    177218                    reltimes.append( None ) 
     219            # Get residuals for each entry. 
    178220            resid = self.planeWaveResiduals( reltimes, relx, rely ) 
    179             print "dbg: names", names 
    180             print "dbg: resid", resid 
    181             print "dbg: modif", self.modifdict 
     221            #print "dbg: names", names 
     222            #print "dbg: resid", resid 
     223            #print "dbg: modif", self.modifdict 
     224            # Compute quality (norm of resid vector, excluding Nones) 
    182225            try: 
    183226                vecnorm = np.sqrt( (resid*resid).sum() ) 
     
    189232                        continue 
    190233                    sqsum += r*r 
    191                 vecnomr = np.sqrt( sqsum ) 
     234                vecnorm = np.sqrt( sqsum ) 
     235            # Store config with minimum residuals. 
    192236            if minvecnorm == None or vecnorm < minvecnorm: 
    193237                minvecnorm = vecnorm 
    194238                optcorrinfo = deepcopy( corrinfo ) 
    195                 print "dbg: copied corrinfo at", vecnorm 
     239                #print "dbg: copied corrinfo at", vecnorm 
     240            # Evaluate improvement to last step. If result is much worse than 
     241            # before, delete this maximum. 
    196242            if lastvecnorm == None: 
    197243                improvement = 0. 
     
    199245                improvement = lastvecnorm - vecnorm 
    200246            if improvement < -5. and modsta != None: 
    201                 # much worse, throw away this corr 
     247                # Much worse, throw away this corr. 
    202248                if len(corrinfo[modsta]) > 1: 
    203249                    corrinfo[modsta] = corrinfo[modsta][1:] 
    204                     print "dbg: deleted entry at", modsta 
     250                    #print "dbg: deleted entry at", modsta 
    205251            lastvecnorm = vecnorm 
    206             print "dbg: resid norm", vecnorm, improvement 
    207             # put first correlation max to last position 
    208             modsta = self.modifyCorrinfo( corrinfo, resid, names ) 
    209         return optcorrinfo 
     252            #print "dbg: resid norm", vecnorm, improvement 
     253            # Put first correlation max to last position and reiterate. 
     254            modsta = self.modifyCorrinfo( corrinfo, resid, names, maxtol ) 
     255            #print "dbg: findOptShift counter", i 
     256            if modsta == None: 
     257                nonecnt += 1 
     258            else: 
     259                nonecnt = 0 
     260            if nonecnt > 3: 
     261                break 
     262        #print "dbg: minimum norm", minvecnorm 
     263        return (optcorrinfo,minvecnorm) 
    210264     
    211265    def planeWaveResiduals( self, reltimes, relx, rely ): 
     266        "Compute residuals of given shift times to a plane wave." 
     267        # Bit complicated due to possible None values. 
    212268        if None in reltimes: 
    213269            rreltimes = [] 
     
    226282        az, azerr, slo, sloerr = computeSlownessFromPicks( np.array(rreltimes), 
    227283            np.array(rrelx), np.array(rrely) ) 
     284        # Now we have slowness and azimuth, let's compute residuals. 
    228285        reltimes_mean = sum(rreltimes)/float(len(rreltimes)) 
    229         # compute beam shifts (as in command 'beam' 
     286        # compute beam shifts (as in command 'beam') 
    230287        radaz = (90.0 - az) / (180./np.pi) 
    231288        kmslo = slo / DEG_TO_KM 
     
    237294        bshifts_mean = sum(bshifts)/float(len(bshifts)) 
    238295        resids = [] 
     296        # Residuals consists of sum (difference) of correlation maxima and beam. 
    239297        for i,relt in enumerate(reltimes): 
    240298            if relt == None: 
     
    245303        return np.array(resids) 
    246304 
    247     def modifyCorrinfo( self, corrinfo, resid, names, maxtol=2. ): 
    248         for i in range(50): 
     305    def modifyCorrinfo( self, corrinfo, resid, names, maxtol=2., loopcnt=50 ): 
     306        """ 
     307        Modify corrinfo for the next iteration. Returns name of changed station. 
     308        resid and names contain the residuals and station names (in the same 
     309        order). 
     310        """ 
     311        for i in range(loopcnt): 
     312            # Determine maximum residual, corresponding station 
    249313            amax = resid.argmax() 
    250314            maxresid = resid[amax] 
    251315            maxsta = names[amax] 
     316            #print "dbg: argmax", amax, maxsta, len(corrinfo[maxsta]), repr(maxresid) 
     317            # Globally count modifications at this station. 
    252318            if not maxsta in self.modifdict: 
    253319                self.modifdict[maxsta] = 0 
    254             # delete entry if max counter reach ed and still over tol 
    255             if self.modifdict[maxsta] >= self.maxmodif \ 
    256                 and len(corrinfo[maxsta]) > 0 and maxresid > maxtol: 
     320            # Delete current (=first) entry if max counter reached or only one 
     321            # entry left and still over tolerance. 
     322            if (self.modifdict[maxsta] >= self.maxmodif \ 
     323                and len(corrinfo[maxsta]) > 0) or len(corrinfo[maxsta]) == 1 \ 
     324                and maxresid > maxtol: 
    257325                if len(corrinfo[maxsta]) > 1: 
     326                    # More than one, reduce by one. 
    258327                    corrinfo[maxsta] = corrinfo[maxsta][1:] 
     328                    #print "dbg: delete entry at", maxsta 
    259329                else: 
     330                    # Delete the last one. 
    260331                    corrinfo[maxsta] = [] 
    261                     print "dbg: *** empty", maxsta 
     332                    #print "dbg: *** empty", maxsta 
    262333                self.modifdict[maxsta] += 1 
    263334                return maxsta 
     335            # Change the station with a large residual below maximum count or 
     336            # a station with a smaller residual and comparable correlation 
     337            # maximum. 
    264338            if self.modifdict[maxsta] < self.maxmodif \ 
    265339                and len(corrinfo[maxsta]) > 1 \ 
    266                 and (corrinfo[maxsta][0][0] > maxtol \ 
     340                and (maxresid > maxtol \ 
    267341                or corrinfo[maxsta][0][1] > 0.95): 
    268                 print "dbg: rotating", maxsta, corrinfo[maxsta][0] 
     342                #print "dbg: rotating", maxsta, corrinfo[maxsta][0] 
    269343                corrinfo[maxsta] = corrinfo[maxsta][1:] + [corrinfo[maxsta][0]] 
    270344                self.modifdict[maxsta] += 1 
    271345                return maxsta 
    272346            resid[amax] = resid[resid.argmin()] 
    273         print "dbg: failed to modify corrinfo" 
     347            #print "dbg: modifyCorrinfo counter", i 
     348        #print "dbg: failed to modify corrinfo" 
    274349        return None 
    275350 
  • SHX/trunk/SeismicHandler/commands/phase.py

    r1158 r1161  
    183183        elif subcmd == 'slowness': 
    184184            pdict = getPhaseSlownesses() 
    185             print "dbg:", pdict 
    186185            pname = self.parameters[1] 
    187186            if pname not in pdict.keys(): 
Note: See TracChangeset for help on using the changeset viewer.