Changeset 1206 for SHX


Ignore:
Timestamp:
19.12.2016 14:38:41 (3 years ago)
Author:
klaus
Message:

ml estimator

File:
1 edited

Legend:

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

    r1204 r1206  
    77import os 
    88import math 
     9import numpy as np 
    910from SeismicHandler.basics.command import BaseCommand 
    1011from SeismicHandler.basics.analysispar import PhaseList, \ 
     
    3233    magnitude determine <magtype> <trclist> <phaselist> <span> 
    3334    magnitude mean <magtype> [<outvar>] 
     35    magnitude estimator <magtype> [<outvar-ml>] [<outvar-mlcorr>] [<picfile>] 
    3436    magnitude clear <magtype> [<station>] 
    3537    magnitude clear_all 
     
    315317            else: 
    316318                set_infodisplay( "determined %d %s values" % (magcnt,magtype) ) 
     319        elif subcmd == 'estimator': 
     320            magtype = self.getMagtype( self.parameters[1] ) 
     321            if magtype != 'ml': 
     322                print "Estimator only for ml" 
     323                return 
     324            trclist = traces_from_list( self.parameters[2] ) 
     325            reslist = [] 
     326            for trc in trclist: 
     327                # station name, skip if already in result list 
     328                sname = "%s.%s" % (trc.stats.network,trc.stats.station) 
     329                if sname not in reslist: 
     330                    degdist = getDistance( trc, epilat, epilon ) 
     331                    mag = maglist.getStationMagnitude(station,magtype) 
     332                    if degdist != None and mag != None: 
     333                        reslist.append( (sname,degdist,mag) ) 
     334            me = MlEstimator() 
     335            magest, corrmag = me.estimateMl( reslist ) 
     336            magest, corrmag = me.estimateMl( reslist, magnstart=magest ) 
     337            if len(self.parameters) > 3: 
     338                magsymb = self.parameters[3][1:] 
     339                self.symbols.set( magsymb, magest ) 
     340            else: 
     341                print "magn", magest 
     342            if len(self.parameters) > 4: 
     343                magsymb = self.parameters[4][1:] 
     344                self.symbols.set( magsymb, corrmag ) 
     345            else: 
     346                print "corr magn", corrmag 
     347            if len(self.parameters) > 5: 
     348                picfile = self.parameters[4].lower() 
     349                me.plot( picfile ) 
    317350        else: 
    318351            raise ShxError( 
     
    472505    return magn 
    473506 
     507class MlEstimator: 
     508 
     509    """Estimate ml magnitude of a preferably large set of station magnitudes. 
     510    Uses output file 'ml_statistics.at' of program ml_statistics.py. 
     511    """ 
     512     
     513    def __init__( self ): 
     514        self.initialized = False 
     515        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 ) 
     520        else: 
     521            self.initialized = False 
     522     
     523    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
     524 
     525    def parseFile( self, fname ): 
     526        self.statcorr = {} 
     527        for line in file(fname): 
     528            if line.startswith("#"): 
     529                continue 
     530            key = "constcorr:" 
     531            if line.startswith(key): 
     532                try: 
     533                    constcorr = float( line[len(key):] ) 
     534                except: 
     535                    return False 
     536                if constcorr != 1: 
     537                    print "MlEstimator: need constant station corrections" 
     538                    return False 
     539                continue 
     540            key = "sensitivity range zeroval:" 
     541            if line.startswith(key): 
     542                try: 
     543                    self.rangestart = float( line[len(key):] ) 
     544                except: 
     545                    print "Cannot read '%s'" % key 
     546                    return False 
     547                continue 
     548            key = "sensitivity range inc:" 
     549            if line.startswith(key): 
     550                try: 
     551                    self.rangeinc = float( line[len(key):] ) 
     552                except: 
     553                    print "Cannot read '%s'" % key 
     554                    return False 
     555                continue 
     556            key = "2nd stage zeroval:" 
     557            if line.startswith(key): 
     558                try: 
     559                    self.secondordercorr_zero = float( line[len(key):] ) 
     560                except: 
     561                    print "Cannot read '%s'" % key 
     562                    return False 
     563                continue 
     564            key = "2nd stage inc:" 
     565            if line.startswith(key): 
     566                try: 
     567                    self.secondordercorr_inc = float( line[len(key):] ) 
     568                except: 
     569                    print "Cannot read '%s'" % key 
     570                    return False 
     571                continue 
     572            if ':' in line: 
     573                try: 
     574                    sstr, cstr = line.split(':') 
     575                    self.statcorr[sstr.strip()] = float( cstr.strip() ) 
     576                except: 
     577                    print "parse error in line:", line 
     578                    return False 
     579        return True 
     580     
     581    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
     582 
     583    def estimateMl( self, values, magnstart=None ): 
     584        """Needs list of triples (station,distance,ml), station like 'GR.BFO', 
     585        distance in deg, ml is uncorrected station magnitude. 
     586        magnstart is a first approximation of the magnitude. If not specified 
     587        derived from the closest stations. 
     588        """ 
     589        if magnstart == None: 
     590            magnstart = self.magnApprox( values ) 
     591            if magnstart == None: 
     592                return (None,None) 
     593        mindist = 30. 
     594        maxdist = self.rangestart + (magnstart-1.)*self.rangeinc 
     595        dv = [] 
     596        mv = [] 
     597        for sta,dist,magn in values: 
     598            if dist < mindist or dist > maxdist: 
     599                continue 
     600            dv.append( dist ) 
     601            if sta in self.statcorr.keys(): 
     602                mv.append( magn-self.statcorr[sta] ) 
     603            else: 
     604                mv.append( magn ) 
     605        dmperkm, zeromag = np.polyfit( np.array(dv), np.array(mv), 1 ) 
     606        corrmag = zeromag-self.secondOrderCorrection(zeromag,dmperkm) 
     607        self.plotdata = (values,mindist,maxdist,zeromag,dmperkm,corrmag) 
     608        return (zeromag,corrmag) 
     609             
     610    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
     611 
     612    def plot( self, fname ): 
     613        if self.plotdata == None: 
     614            return 
     615        values,mindist,maxdist,zeromag,dmperkm,corrmag = self.plotdata 
     616        fig = plt.figure( facecolor='white' ) 
     617        ax = fig.add_subplot(111) 
     618        ax.set_xlabel( "distance (km)" ) 
     619        ax.set_ylabel( "Magnitude" ) 
     620        ax.set_title( "ml: %4.2f  corr ml: %4.2f" % (zeromag,corrmag) ) 
     621        ax.set_xlim( 0., 500. ) 
     622        ax.set_ylim( 0.4, 4.5 ) 
     623        labels, dists, mags = zip(*values) 
     624        maxd = max( dists ) 
     625        ax.add_patch( 
     626            patches.Rectangle( 
     627                (0.,0.),mindist,5.,facecolor='yellow', 
     628                edgecolor='white',alpha=0.5 
     629            ) 
     630        ) 
     631        ax.add_patch( 
     632            patches.Rectangle( 
     633                (maxdist,0.),500.-drange[1],5.,facecolor='yellow', 
     634                edgecolor='white',alpha=0.5 
     635            ) 
     636        ) 
     637        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): 
     642            plt.annotate( lab, xy=(x,y), fontsize=9 ) 
     643        plt.draw() 
     644        plt.savefig( fname ) 
     645        plt.close( fig ) 
     646         
     647 
     648    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
     649 
     650    def magnApprox( self, val ): 
     651        "Estimate magnitude from closest 10% of stations." 
     652        if len(val) < 20: 
     653            return None 
     654        sta, dist, mag = zip(*val) 
     655        limdist = sorted(dist)[len(dist)/10] 
     656        mags = [] 
     657        for sta,dist,mag in val: 
     658            if dist <= limdist: 
     659                mags.append( mag ) 
     660        if len(mags) > 4: 
     661            mags.remove( min(mags) ) 
     662            mags.remove( max(mags) ) 
     663        return (sum(mags)/float(len(mags))) 
     664 
     665    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
     666 
     667    def secondOrderCorrection( self, magn, dmperkm ): 
     668        "Returns 2nd order correction for magnitude." 
     669        if self.secondordercorr_zero == None: 
     670            return 0. 
     671        else: 
     672            return (self.secondordercorr_zero \ 
     673                + self.secondordercorr_inc*dmperkm) 
     674 
Note: See TracChangeset for help on using the changeset viewer.