Changeset 1150


Ignore:
Timestamp:
30.01.2016 18:42:29 (4 years ago)
Author:
klaus
Message:

phase plot; reset entry; bugfixes

Location:
SHX/trunk/SeismicHandler
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed
  • SHX/trunk/SeismicHandler/basics/analysispar.py

    r1149 r1150  
    269269        for station in self.magndict.keys(): 
    270270            self.clearStationMagnitude( station, magtype ) 
     271     
     272    def clearAll( self ): 
     273        self.magndict = {} 
    271274     
    272275    def getStations( self ): 
  • SHX/trunk/SeismicHandler/commands/magnitude.py

    r1149 r1150  
    1414from SeismicHandler.basics.error import ShxError 
    1515from SeismicHandler.basics import timeit 
    16 from SeismicHandler.modules.traces import traces_from_list 
     16from SeismicHandler.modules.traces import traces_from_list, getDistance 
    1717from obspy.sh.core import toUTCDateTime, fromUTCDateTime 
    18 from obspy.core.util.geodetics import gps2DistAzimuth 
    1918import matplotlib.pyplot as plt 
    2019 
     
    3130    magnitude mean <magtype> [<outvar>] 
    3231    magnitude clear <magtype> [<station>] 
     32    magnitude clear_all 
    3333    magnitude set <magtype> <station> <value> 
    3434    magnitude plot <magtype> 
     
    6868    legal_subcmds = { 
    6969        'dump'      : 1, 
     70        'clear_all' : 1, 
    7071        'plot'      : 2, 
    7172        'determine' : 5, 
     
    104105            else: 
    105106                print mm 
     107        elif subcmd == 'clear_all': 
     108            maglist.clearAll() 
    106109        elif subcmd == 'clear': 
    107110            magtype = self.getMagtype( self.parameters[1] ) 
     
    187190                    spans.append( (stime,etime) ) 
    188191                # determine maximum magnitude in all spans 
    189                 degdist = self.getDistance( trc, epilat, epilon ) 
     192                degdist = getDistance( trc, epilat, epilon ) 
    190193                if degdist == None: 
    191194                    print "No distance for '%s'" % sname 
     
    231234                reltime = lphase.picktime - trc.stats.starttime \ 
    232235                    - trc.get_info('t-origin') 
     236                if reltime == None: 
     237                    print "dbg: *** no reltime on", trc 
     238                    return 
    233239                # manual phase -> window after pick 
    234240                return (reltime,(reltime+wdwlen)) 
     
    243249        return (None,None) 
    244250     
    245     def getDistance( self, trc, lat, lon ): 
    246         stations = Stations() 
    247         sname = "%s.%s" % (trc.stats.network,trc.stats.station) 
    248         metatime = trc.stats.starttime \ 
    249             + (trc.stats.endtime-trc.stats.starttime)/2. 
    250         try: 
    251             r = stations[(sname,metatime)] 
    252             slat = float( r.latitude ) 
    253             slon = float( r.longitude ) 
    254         except: 
    255             return None 
    256         dist, az, baz = gps2DistAzimuth( lat, lon, slat, slon ) 
    257         return (dist/1000./DEG_TO_KM) 
    258      
    259251    def determineMagnitude( self, trc, magtype, stime, etime, degdist ): 
    260252        if magtype == 'ml': 
     
    286278        for t in hlist: 
    287279            sl = trc.slice_relative( stime, etime ) 
    288             ampl.append( max( abs(max(sl.data)), abs(min(sl.data)) ) ) 
     280            #ampl.append( max( abs(max(sl.data)), abs(min(sl.data)) ) ) 
     281            mean = sl.data.mean() 
     282            ampl.append( max(sl.data.max()-mean, mean-sl.data.min()) ) 
    289283        ampl = math.sqrt( ampl[0]*ampl[0] + ampl[1]*ampl[1] ) 
    290284        return magn_ml_analytic( ampl, degdist*DEG_TO_KM ) 
     
    302296            if station in labels: 
    303297                continue 
    304             dist = self.getDistance( trc, epilat, epilon ) 
     298            dist = getDistance( trc, epilat, epilon ) 
    305299            if dist == None: 
    306300                continue 
  • SHX/trunk/SeismicHandler/commands/param.py

    r1139 r1150  
    9999            ui_event( "updateparams" ) 
    100100        elif subcmd == 'reset': 
    101             analysispar.reset() 
     101            analysispar.resetParams() 
    102102            ui_event( "updateparams" ) 
    103103        elif subcmd == 'savetraces': 
  • SHX/trunk/SeismicHandler/commands/phase.py

    r1146 r1150  
    55#    http://www.seismic-handler.org/portal/wiki/Shx/LicenseTerms 
    66 
     7import os 
    78from SeismicHandler.basics.command import BaseCommand 
    8 from SeismicHandler.basics.analysispar import PhaseList 
     9from SeismicHandler.basics.analysispar import PhaseList, AnalysisPar 
    910from SeismicHandler.basics.error import ShxError 
    1011from SeismicHandler.basics import timeit 
    11 from SeismicHandler.modules.traces import addNetcode, number_of_visible_traces 
     12from SeismicHandler.modules.traces import addNetcode, number_of_visible_traces,\ 
     13    getDistance 
    1214from obspy.sh.core import toUTCDateTime, fromUTCDateTime 
     15from SeismicHandler.config import Settings 
     16import matplotlib.pyplot as plt 
     17import numpy as np 
     18 
     19DEG_TO_KM = 111.19 
    1320 
    1421provides = {"phase": "phase"} 
     
    1623    """ 
    1724    URI:http://www.seismic-handler.org/portal/wiki/ShPhase 
     25     
     26    PHASE DEFINE <station> <picktime> [<phasename>] [<picktype>] [<comp>] 
     27    PHASE LIST <station> [<picktype>] |<comp>] [<minweight>] 
     28    PHASE DEFAULT_PHASE <phase> 
     29    PHASE STATIONS 
     30    PHASE CLEAR [<phase>] [<station>] [<picktype>] 
     31    PHASE CLEAR_ALL 
     32    PHASE DUMP 
     33    PHASE PLOT [<outfile>] 
    1834    """ 
    1935    numberOfParameters = [1,2,3,4,5,6] 
     
    148164                        fromUTCDateTime(phase.picktime), 
    149165                        phaselist.picktypeName(phase.picktype),phase.comp) 
     166        elif subcmd == 'plot': 
     167            if len(self.parameters) > 1: 
     168                outfile = self.parameters[1] 
     169            else: 
     170                outfile = self.tmpPlotfile() 
     171            self.createPlot( outfile ) 
     172            try: 
     173                dspprog = Settings.config.misc.display_prog 
     174            except: 
     175                dspprog = 'display' 
     176            os.system( "%s %s &" % (dspprog,outfile) ) 
    150177        else: 
    151178            raise ShxError( "Command 'phase', unknown subcommand '%s'" % subcmd, 
    152179                status=1111 ) 
     180 
     181    def tmpPlotfile( self ): 
     182        return os.path.join( Settings.config.paths.scratch[0], 
     183            'tmpphasefig_%d.png' % os.getpid() ) 
     184     
     185    def createPlot( self, outfile ): 
     186        ap = AnalysisPar() 
     187        epilat = ap.getValue( 'epi_latitude' ) 
     188        epilon = ap.getValue( 'epi_longitude' ) 
     189        origin = ap.getValue( 'origin_time' ) 
     190        title = "Distance vs Time" 
     191        fig = plt.figure( facecolor='white' ) 
     192        ax = fig.add_subplot(111) 
     193        ax.set_title( title ) 
     194        (times, dists, labels, colors, phases), lineinfo = \ 
     195            self.getPhasePlotData( origin, epilat, epilon ) 
     196        islocal = 'Pg' in phases or 'Pn' in phases or 'Sg' in phases 
     197        if islocal: 
     198            dists = np.array(dists)*DEG_TO_KM 
     199            ax.set_ylabel( "distance (km)" ) 
     200        else: 
     201            ax.set_ylabel( "distance (deg)" ) 
     202        ax.set_xlabel( "time since origin (s)" ) 
     203        mind = min(dists) 
     204        maxd = max(dists)*1.1 
     205        for pname,slope,color in lineinfo: 
     206            if islocal: 
     207                slope /= DEG_TO_KM 
     208                lab = "%s (%4.2f km/s)" % (pname,(1./slope)) 
     209            else: 
     210                lab = "%s (%g s/deg)" % (pname,slope) 
     211            ax.plot( [mind*slope, maxd*slope], [mind, maxd], c=color ) 
     212            plt.annotate( lab, xy=(maxd*slope,maxd*1.01), fontsize=9, 
     213                horizontalalignment='center', color=color ) 
     214        ax.scatter( times, dists, c=colors, edgecolors=colors ) 
     215        for lab,x,y in zip(labels,times,dists): 
     216            plt.annotate( lab, xy=(x,y), fontsize=9 ) 
     217        plt.draw() 
     218        plt.savefig( outfile, facecolor=fig.get_facecolor() ) 
     219 
     220    def getPhasePlotData( self, origin, lat, lon ): 
     221        colorlist = [(1,0,0),(0,0,1),(0,1,0),(1,1,0),(0,1,1)] 
     222        colidx = 0 
     223        coldict = {} 
     224        phasedict = {}  # collect dist,tt per phase type 
     225        phaselist = PhaseList() 
     226        tinfo = [] # separate list for theo phases, so they can be plotted first 
     227        info = [] 
     228        manual = phaselist.translatePicktype( 'manual' ) 
     229        auto = phaselist.translatePicktype( 'auto' ) 
     230        theo = phaselist.translatePicktype( 'theo' ) 
     231        for station in phaselist.getStations(): 
     232            haspick = False 
     233            linfo = [] 
     234            ltinfo = [] 
     235            dist = getDistance( station, lat, lon ) 
     236            if dist == None: 
     237                continue 
     238            for phase in phaselist.getPhaseList(station): 
     239                if phase.picktype not in (manual,auto,theo): 
     240                    continue 
     241                if phase.picktype != theo: 
     242                    haspick = True 
     243                ptime = phase.picktime - origin 
     244                if phase.picktype == theo: 
     245                    lab = '' 
     246                    col = (0.78,0.78,0.78) 
     247                    ltinfo.append( (ptime,dist,lab,col,phase.name) ) 
     248                else: 
     249                    try: 
     250                        lab = station.split('.')[1] 
     251                    except: 
     252                        lab = station 
     253                    if not phase.name in coldict: 
     254                        if colidx == len(colorlist): 
     255                            colidx = 0 
     256                        coldict[phase.name] = colorlist[colidx] 
     257                        colidx += 1 
     258                    col = coldict[phase.name] 
     259                    if not phase.name in phasedict.keys(): 
     260                        phasedict[phase.name] = [] 
     261                    phasedict[phase.name].append( (dist,ptime) ) 
     262                    linfo.append( (ptime,dist,lab,col,phase.name) ) 
     263            if haspick: 
     264                info += linfo 
     265                tinfo += ltinfo 
     266        # slope and line fit for each phase 
     267        lineinfo = [] 
     268        for pname in phasedict.keys(): 
     269            slope = 0. 
     270            for dist,ptime in phasedict[pname]: 
     271                slope += ptime/dist 
     272            slope /= float(len(phasedict[pname])) 
     273            lineinfo.append( (pname,slope,coldict[pname]) ) 
     274        return (zip(*(tinfo+info)),lineinfo) 
     275                 
     276             
  • SHX/trunk/SeismicHandler/commands/sort_by_distance.py

    r1083 r1150  
    77from SeismicHandler.basics.command import BaseCommand 
    88from SeismicHandler.core import Traces 
    9 from SeismicHandler.modules.traces import Traces as BaseTraces, traces_from_list 
     9from SeismicHandler.modules.traces import Traces as BaseTraces, \ 
     10    traces_from_list, getDistance 
    1011from SeismicHandler.basics.error import ShxError 
    1112from SeismicHandler.basics import timeit 
    1213from SeismicHandler.modules.stations import Stations 
    13 from obspy.core.util.geodetics import gps2DistAzimuth 
    1414 
    1515 
     
    4646        stations = Stations() 
    4747        for t in lst: 
    48             sname = "%s.%s.%s.%s" % (t.stats.network, t.stats.station, 
    49                 t.stats.location, t.stats.channel) 
    50             metatime = t.stats.starttime + (t.stats.endtime-t.stats.starttime)/2 
    51             try: 
    52                 r = stations[(sname, metatime)] 
    53             except: 
    54                 r = None 
    55             try: 
    56                 slat = float( r.latitude ) 
    57                 slon = float( r.longitude ) 
    58             except: 
     48            dist = getDistance( t, lat, lon ) 
     49            if dist == None: 
    5950                continue 
    60             dist, az, baz = gps2DistAzimuth( slat, slon, lat, lon ) 
    6151            t.shx._parent.remove(t) 
    6252            affected.append( (t,dist) ) 
  • SHX/trunk/SeismicHandler/commands/theophase.py

    r1148 r1150  
    1010from SeismicHandler.basics import timeit 
    1111from SeismicHandler.basics.error import ShxError 
    12 from SeismicHandler.modules.traces import traces_from_list 
     12from SeismicHandler.modules.traces import traces_from_list, getDistance 
    1313from SeismicHandler.modules.traveltime import LocalTravelTime, TaupiTravelTime 
    1414from SeismicHandler.basics.analysispar import AnalysisPar, PhaseList 
    1515from SeismicHandler.modules.stations import Stations 
    16 from obspy.core.util.geodetics import gps2DistAzimuth 
    1716from obspy.sh.core import toUTCDateTime 
    1817 
     
    132131        distdict = {} 
    133132        for trc in traces: 
     133            dist = getDistance( trc, epilat, epilon ) 
    134134            sname = "%s.%s" % (trc.stats.network,trc.stats.station) 
    135             if sname in distdict.keys(): 
    136                 continue 
    137             metatime = trc.stats.starttime \ 
    138                 + (trc.stats.endtime-trc.stats.starttime)/2 
    139             try: 
    140                 r = stations[(sname,metatime)] 
    141                 slat = float( r.latitude ) 
    142                 slon = float( r.longitude ) 
    143             except: 
    144                 continue 
    145             dist, az, baz = gps2DistAzimuth( epilat, epilon, slat, slon ) 
    146135            if distunit == 'km': 
    147                 dist /= 1000. 
    148             else: 
    149                 dist /= (1000.*DEG_TO_KM) 
     136                dist *= DEG_TO_KM 
    150137            distdict[sname] = dist 
    151138 
  • SHX/trunk/SeismicHandler/modules/traces.py

    r1149 r1150  
    1313from SeismicHandler.config import Settings, set_dsptrcs, set_tottrcs 
    1414from SeismicHandler.modules.stations import Stations 
     15from obspy.core.util.geodetics import gps2DistAzimuth 
    1516 
    1617 
     
    1920META_STATUS_COMPLETE = 2 
    2021 
     22DEG_TO_KM = 111.19 
    2123 
    2224class Traces(object): 
     
    449451    return "%s.%s" % (netcode,upstat) 
    450452 
     453def getDistance( trc, lat, lon ): 
     454    """ 
     455    Compute distance in degrees from epicenter to station of trace. 
     456    'trc' is an obspy trace or a station name (at least NET.STATION). 
     457    """ 
     458    if trc == None: 
     459        return None 
     460    if hasattr(trc,'lower'): 
     461        # trc is a string -> station name 
     462        sname = trc 
     463        # have to take time from first trace on display 
     464        try: 
     465            trc = Traces().traces.visible[0] 
     466        except: 
     467            return None 
     468    else: 
     469        # trc is a obspy trace, take station name from trace metadata 
     470        sname = "%s.%s" % (trc.stats.network,trc.stats.station) 
     471    metatime = trc.stats.starttime \ 
     472        + (trc.stats.endtime-trc.stats.starttime)/2. 
     473    stations = Stations() 
     474    try: 
     475        r = stations[(sname,metatime)] 
     476        slat = float( r.latitude ) 
     477        slon = float( r.longitude ) 
     478    except: 
     479        return None 
     480    dist, az, baz = gps2DistAzimuth( lat, lon, slat, slon ) 
     481    return (dist/1000./DEG_TO_KM) 
     482     
    451483 
    452484if __name__ == "__main__": 
  • SHX/trunk/SeismicHandler/modules/wx_.py

    r1149 r1150  
    321321        if trc == None: 
    322322            self.clearWindow() 
     323            return 
    323324        if self._s3_showalltraces: 
    324325            s3_recent = self._s3_origtrcname 
     
    17881789        marginfactor = 0.3 
    17891790        mf2 = (1 - marginfactor) * 2 
    1790         th2 = self.traceheight / 2 
     1791        th2 = self.traceheight / 2 + 1 
    17911792        offset = midpoint - th2 - self.Scrolled 
    17921793        if self.traceOrder in [0, 1]: 
    17931794            self._drawXor( 'rect', (start[0], offset + th2 * marginfactor, 
    1794                 end[0] - start[0], th2 * mf2), dc ) 
     1795                end[0] - start[0], th2 * mf2 + 1), dc ) 
    17951796        elif self.traceOrder in [2, 3]: 
    17961797            self._drawXor( 'rect', (offset + th2 * marginfactor, start[1], 
     
    21002101        self.addEntry( phaseMenu, 'Delete all Theo Phases', 
    21012102            'Delete all theoretical phases', self.OnDeleteTheoPhase ) 
     2103        phaseMenu.AppendSeparator() 
     2104        self.addEntry( phaseMenu, 'Plot Phases', 
     2105            'Create distance vs travel time plot', self.OnPhasePlot ) 
    21022106        # magnitude menu 
    21032107        magnMenu = wx.Menu() 
     
    21302134            self.OnRecoverTracesAndParams ) 
    21312135        controlMenu.AppendSeparator() 
     2136        self.addEntry( controlMenu, 'Clear All Parameters and Traces', 
     2137            'Deletes traces, phase, magnitudes and resets parameters.', 
     2138            self.OnClearAllParameters ) 
     2139        controlMenu.AppendSeparator() 
    21322140        self.addEntry( controlMenu, 'Open Parameter Window', 
    21332141            'Open parameter dialog', self.OnOpenParams ) 
     
    25422550        "plotterWindow: control menu entry." 
    25432551        _sendShCommand( 'TRACES_PARAMS_RECOVER' ) 
    2544      
     2552         
     2553    def OnClearAllParameters( self, e ): 
     2554        _sendShCommand( 'shx_menu_reset_all' ) 
     2555         
    25452556    def OnRecoverTracesAndParams( self, e ): 
    25462557        "plotterWindow: control menu entry." 
     
    26232634        "plotterWindow: phase menu entry." 
    26242635        _sendShCommand( "phase clear ;;; theo\ntheophase iasp91 all" ) 
     2636     
     2637    def OnPhasePlot( self, e ): 
     2638        _sendShCommand( "phase plot" ) 
    26252639 
    26262640    def OnDeleteTheoPhase( self, e ): 
Note: See TracChangeset for help on using the changeset viewer.