Changeset 1207 for SHX


Ignore:
Timestamp:
19.12.2016 17:16:17 (3 years ago)
Author:
klaus
Message:

bugfixes and ml estimator script

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

Legend:

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

    r1206 r1207  
    1818from SeismicHandler.modules.traces import traces_from_list, getDistance 
    1919import matplotlib.pyplot as plt 
     20import matplotlib.patches as patches 
    2021from obspy.io.sh.core import to_utcdatetime as toUTCDateTime 
    2122from obspy.io.sh.core import from_utcdatetime as fromUTCDateTime 
     
    109110     Delete ml magnitude on station gr.gra1. 
    110111    """ 
    111     numberOfParameters = [1,2,3,4,5] 
     112    numberOfParameters = [1,2,3,4,5,6] 
    112113    parameterQueries = [ 
    113114        { 
     
    150151        'clear'     : 2, 
    151152        'set'       : 4, 
     153        'estimator' : 3, 
    152154    } 
    153155     
     
    318320                set_infodisplay( "determined %d %s values" % (magcnt,magtype) ) 
    319321        elif subcmd == 'estimator': 
     322            ap = AnalysisPar() 
     323            epilat = ap.getValue( 'epi_latitude' ) 
     324            epilon = ap.getValue( 'epi_longitude' ) 
     325            if epilat == None or epilon == None: 
     326                print "No epicenter -> no magnitude" 
     327                return 
    320328            magtype = self.getMagtype( self.parameters[1] ) 
    321329            if magtype != 'ml': 
     
    324332            trclist = traces_from_list( self.parameters[2] ) 
    325333            reslist = [] 
     334            stalist = [] 
    326335            for trc in trclist: 
    327336                # station name, skip if already in result list 
    328337                sname = "%s.%s" % (trc.stats.network,trc.stats.station) 
    329                 if sname not in reslist: 
     338                if sname not in stalist: 
    330339                    degdist = getDistance( trc, epilat, epilon ) 
    331                     mag = maglist.getStationMagnitude(station,magtype) 
     340                    mag = maglist.getStationMagnitude(sname,magtype) 
    332341                    if degdist != None and mag != None: 
    333                         reslist.append( (sname,degdist,mag) ) 
     342                        reslist.append( (sname,degdist*DEG_TO_KM,mag) ) 
     343                        stalist.append( sname ) 
    334344            me = MlEstimator() 
    335345            magest, corrmag = me.estimateMl( reslist ) 
     
    346356                print "corr magn", corrmag 
    347357            if len(self.parameters) > 5: 
    348                 picfile = self.parameters[4].lower() 
     358                picfile = self.parameters[5].lower() 
    349359                me.plot( picfile ) 
    350360        else: 
     
    514524        self.initialized = False 
    515525        self.plotdata = None 
    516         filename = os.path.join( Settings.config.paths.private, 
    517             'ml_statistics.dat' ) 
    518         if os.path.exists(filename): 
    519             self.initialized = self.parseFile( filename ) 
     526        for pname in Settings.config.paths.private: 
     527            filename = os.path.join( pname, 'ml_statistics.dat' ) 
     528            if os.path.exists(filename): 
     529                self.initialized = self.parseFile( filename ) 
     530                break 
    520531        else: 
    521532            self.initialized = False 
     
    621632        ax.set_xlim( 0., 500. ) 
    622633        ax.set_ylim( 0.4, 4.5 ) 
    623         labels, dists, mags = zip(*values) 
     634        labels, dists, umags = zip(*values) 
     635        # apply corrections umags -> cmags 
     636        cmags = [] 
     637        for sta,mag in zip(labels,umags): 
     638            if sta in self.statcorr.keys(): 
     639                cmags.append( mag-self.statcorr[sta] ) 
     640            else: 
     641                cmags.append( mag ) 
    624642        maxd = max( dists ) 
    625643        ax.add_patch( 
     
    631649        ax.add_patch( 
    632650            patches.Rectangle( 
    633                 (maxdist,0.),500.-drange[1],5.,facecolor='yellow', 
     651                (maxdist,0.),500.-maxdist,5.,facecolor='yellow', 
    634652                edgecolor='white',alpha=0.5 
    635653            ) 
    636654        ) 
    637655        ax.plot( [0.,maxd], [zeromag,(zeromag+dmperkm*maxd)], c='r' ) 
    638         ax.scatter( dists, mags ) 
    639         if min(mags) < 0.: 
    640             ax.set_ylim( 0., 1.1*max(mags) ) 
    641         for lab,x,y in zip(labels,dists,mags): 
     656        ax.scatter( dists, umags, c='#eeeeee', edgecolor='' ) 
     657        ax.scatter( dists, cmags, c='blue' ) 
     658        if min(cmags) < 0.: 
     659            ax.set_ylim( 0., 1.1*max(cmags) ) 
     660        for lab,x,y in zip(labels,dists,cmags): 
    642661            plt.annotate( lab, xy=(x,y), fontsize=9 ) 
    643662        plt.draw() 
Note: See TracChangeset for help on using the changeset viewer.