Changeset 1160 for SHX


Ignore:
Timestamp:
05.02.2016 18:28:04 (4 years ago)
Author:
klaus
Message:

corrpick, evaluating residuals, not finished

Location:
SHX/trunk/SeismicHandler
Files:
2 edited

Legend:

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

    r1158 r1160  
    99""" 
    1010 
    11 # Standard SeismicHandler libraries to import. 
     11import numpy as np 
     12from copy import deepcopy 
    1213from SeismicHandler.basics.command import BaseCommand 
    1314from SeismicHandler.modules.traces import traces_from_list 
     15from SeismicHandler.modules.stations import Stations 
     16from SeismicHandler.modules.seismics import absToRelPositions, \ 
     17    computeSlownessFromPicks 
    1418from SeismicHandler.basics.messages import log_message 
    1519from SeismicHandler.basics.analysispar import PhaseList 
     
    1721from SeismicHandler.basics import timeit 
    1822from obspy.signal.cross_correlation import xcorr 
     23 
     24DEG_TO_KM = 111.19 
    1925 
    2026# Declare which command is implemented here. 
     
    6369        "Correlation picker." 
    6470 
     71        nummax = 3 
     72        self.modifdict = {} 
     73        self.maxmodif = 2*nummax 
     74 
    6575        traces = traces_from_list(self.parameters[0]) 
    6676        if len(traces) != 1: 
     
    7181        _end = float( self.parameters[2] ) 
    7282        _shift = float( self.parameters[3] ) 
    73         print "dbg: parameters", _start, _end, _shift 
    7483        wdwlen = _end - _start 
    7584        trace_sname = "%s.%s" % (trace.stats.network,trace.stats.station) 
     
    9099        shiftlength = int( _shift / trace.stats.delta )/2 - 1 
    91100        auto = self.phaselist.translatePicktype( 'auto' ) 
     101        corrinfo = {} 
    92102        for trc in traces_from_list("all"): 
    93103            sname = "%s.%s" % (trc.stats.network,trc.stats.station) 
     
    96106            if trc.stats.channel != trace.stats.channel: 
    97107                continue 
    98             print "dbg: len", len(trc.data), shiftlength 
    99108            idx, val, cfct = xcorr( corrtrc, trc.slice(shiftstart,shiftend),  
    100109                shiftlength, full_xcorr=True ) 
    101             idx = cfct.argmax() - len(cfct)/2 
    102             print "dbg:", sname, (trace.stats.delta*idx), val 
     110            corrinfo[sname] = [] 
     111            maxcorrval = None 
     112            for m in range(nummax): 
     113                amax = cfct.argmax() 
     114                idx = amax - len(cfct)/2 
     115                corrval = cfct[amax] 
     116                if maxcorrval == None: 
     117                    maxcorrval = corrval 
     118                reltime = -trace.stats.delta*idx 
     119                corrinfo[sname].append( (reltime,(corrval/maxcorrval)) ) 
     120                minval = cfct.min() 
     121                cfct[amax] = minval  # delete this maximum 
     122                try: 
     123                    cfct[amax-1] = minval 
     124                    cfct[amax+1] = minval 
     125                    cfct[amax-2] = minval 
     126                    cfct[amax+2] = minval 
     127                except: 
     128                    pass 
     129        corrinfo = self.findOptShifts( corrinfo ) 
     130        for sname in corrinfo.keys(): 
     131            reltime = corrinfo[sname][0][0] 
    103132            self.phaselist.clearPhasenameStation( phase, sname, auto ) 
    104             reltime = -trace.stats.delta*idx 
     133            if reltime == None: 
     134                continue 
    105135            self.phaselist.addPhasePick( sname, phasetime+reltime, name=phase, 
    106136                picktype=auto ) 
     
    119149            return (phase.name,phase.picktime) 
    120150        return (None,None) 
    121  
     151     
     152    def findOptShifts( self, corrinfo ): 
     153        # find station locations 
     154        optcorrinfo = None 
     155        minvecnorm = None 
     156        lastvecnorm = None 
     157        stations = Stations() 
     158        lats = [] 
     159        lons = [] 
     160        names = corrinfo.keys() 
     161        for station in names: 
     162            try: 
     163                r = stations[(station,None)] 
     164                lat = float( r.latitude ) 
     165                lon = float( r.longitude ) 
     166            except: 
     167                lat = lon = None 
     168            lats.append( lat ) 
     169            lons.append( lon ) 
     170        relx, rely = absToRelPositions( lons, lats ) 
     171        for i in range(50): 
     172            reltimes = [] 
     173            for s in names: 
     174                if len(corrinfo[s]) > 0: 
     175                    reltimes.append( corrinfo[s][0][0] ) 
     176                else: 
     177                    reltimes.append( None ) 
     178            resid = self.planeWaveResiduals( reltimes, relx, rely ) 
     179            print "dbg: names", names 
     180            print "dbg: resid", resid 
     181            print "dbg: modif", self.modifdict 
     182            try: 
     183                vecnorm = np.sqrt( (resid*resid).sum() ) 
     184            except: 
     185                # None's in it 
     186                sqsum = 0. 
     187                for r in resid: 
     188                    if r == None: 
     189                        continue 
     190                    sqsum += r*r 
     191                vecnomr = np.sqrt( sqsum ) 
     192            if minvecnorm == None or vecnorm < minvecnorm: 
     193                minvecnorm = vecnorm 
     194                optcorrinfo = deepcopy( corrinfo ) 
     195                print "dbg: copied corrinfo at", vecnorm 
     196            if lastvecnorm == None: 
     197                improvement = 0. 
     198            else: 
     199                improvement = lastvecnorm - vecnorm 
     200            if improvement < -5. and modsta != None: 
     201                # much worse, throw away this corr 
     202                if len(corrinfo[modsta]) > 1: 
     203                    corrinfo[modsta] = corrinfo[modsta][1:] 
     204                    print "dbg: deleted entry at", modsta 
     205            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 
     210     
     211    def planeWaveResiduals( self, reltimes, relx, rely ): 
     212        if None in reltimes: 
     213            rreltimes = [] 
     214            rrelx = [] 
     215            rrely = [] 
     216            for i,relt in enumerate(reltimes): 
     217                if relt == None: 
     218                    continue 
     219                rreltimes.append( relt ) 
     220                rrelx.append( relx[i] ) 
     221                rrely.append( rely[i] ) 
     222        else: 
     223            rreltimes = reltimes 
     224            rrelx = relx 
     225            rrely = rely 
     226        az, azerr, slo, sloerr = computeSlownessFromPicks( np.array(rreltimes), 
     227            np.array(rrelx), np.array(rrely) ) 
     228        reltimes_mean = sum(rreltimes)/float(len(rreltimes)) 
     229        # compute beam shifts (as in command 'beam' 
     230        radaz = (90.0 - az) / (180./np.pi) 
     231        kmslo = slo / DEG_TO_KM 
     232        bshifts = [] 
     233        for i,relt in enumerate(reltimes): 
     234            bshifts.append( 
     235                (relx[i]*np.cos(radaz) + rely[i]*np.sin(radaz)) * kmslo 
     236            ) 
     237        bshifts_mean = sum(bshifts)/float(len(bshifts)) 
     238        resids = [] 
     239        for i,relt in enumerate(reltimes): 
     240            if relt == None: 
     241                resid = None 
     242            else: 
     243                resid = abs( (bshifts[i]-bshifts_mean) + (relt-reltimes_mean) ) 
     244            resids.append( resid ) 
     245        return np.array(resids) 
     246 
     247    def modifyCorrinfo( self, corrinfo, resid, names, maxtol=2. ): 
     248        for i in range(50): 
     249            amax = resid.argmax() 
     250            maxresid = resid[amax] 
     251            maxsta = names[amax] 
     252            if not maxsta in self.modifdict: 
     253                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: 
     257                if len(corrinfo[maxsta]) > 1: 
     258                    corrinfo[maxsta] = corrinfo[maxsta][1:] 
     259                else: 
     260                    corrinfo[maxsta] = [] 
     261                    print "dbg: *** empty", maxsta 
     262                self.modifdict[maxsta] += 1 
     263                return maxsta 
     264            if self.modifdict[maxsta] < self.maxmodif \ 
     265                and len(corrinfo[maxsta]) > 1 \ 
     266                and (corrinfo[maxsta][0][0] > maxtol \ 
     267                or corrinfo[maxsta][0][1] > 0.95): 
     268                print "dbg: rotating", maxsta, corrinfo[maxsta][0] 
     269                corrinfo[maxsta] = corrinfo[maxsta][1:] + [corrinfo[maxsta][0]] 
     270                self.modifdict[maxsta] += 1 
     271                return maxsta 
     272            resid[amax] = resid[resid.argmin()] 
     273        print "dbg: failed to modify corrinfo" 
     274        return None 
    122275 
    123276# Standard main part (same for all commands). 
  • SHX/trunk/SeismicHandler/modules/wx_.py

    r1159 r1160  
    193193    def OnKeyPress( self, evt ): 
    194194        "magnifyCanvas: key event: emulate all shortcuts of main window." 
    195         char = chr( evt.GetUniChar() ).lower() 
     195        localtrans = { 
     196            314 : '<-', 
     197            315 : '-^', 
     198            316 : '->', 
     199            317 : '-v', 
     200        } 
     201        try: 
     202            char = chr( evt.GetUniChar() ).lower() 
     203        except: 
     204            try: 
     205                char = localtrans[evt.GetUniChar()] 
     206            except: 
     207                print "dbg: unichar not recognized", repr(evt.GetUniChar()) 
     208                return 
    196209        shortcut = { 
    197210            'b': plotter.OnBeam, 
     
    214227            '3': plotter.OnShowAllInMagnify, 
    215228            '8': plotter.OnFilterBP_1_8, 
     229            '-^': self.keyZoomUp, 
     230            '-v': self.keyZoomDown, 
     231            '->': self.keyZoomGrow, 
     232            '<-': self.keyZoomShrink, 
    216233        } 
    217234        shortcut_control = { 
     
    232249            if char in shortcut.keys(): 
    233250                shortcut[char](None) 
     251     
     252    def keyZoomGrow( self, e ): 
     253        plotter.canvas.moveZoomWindow( 'grow' ) 
     254    def keyZoomShrink( self, e ): 
     255        plotter.canvas.moveZoomWindow( 'shrink' ) 
     256    def keyZoomUp( self, e ): 
     257        plotter.canvas.moveZoomWindow( 'up' ) 
     258    def keyZoomDown( self, e ): 
     259        plotter.canvas.moveZoomWindow( 'down' ) 
    234260 
    235261    def refresh( self ): 
Note: See TracChangeset for help on using the changeset viewer.