source: SHX/trunk/SeismicHandler/commands/magnitude.py @ 1207

Revision 1207, 25.7 KB checked in by klaus, 3 years ago (diff)

bugfixes and ml estimator script

Line 
1# -*- coding: utf-8 -*-
2
3#    This file is part of Seismic Handler eXtended (SHX). For terms of use and
4#    license information please see license.txt and visit
5#    http://www.seismic-handler.org/portal/wiki/Shx/LicenseTerms
6
7import os
8import math
9import numpy as np
10from SeismicHandler.basics.command import BaseCommand
11from SeismicHandler.basics.analysispar import PhaseList, \
12    MagnitudeList, AnalysisPar
13from SeismicHandler.config import Settings
14from SeismicHandler.modules.stations import Stations
15from SeismicHandler.basics.error import ShxError
16from SeismicHandler.basics import timeit
17from SeismicHandler.basics.messages import set_infodisplay
18from SeismicHandler.modules.traces import traces_from_list, getDistance
19import matplotlib.pyplot as plt
20import matplotlib.patches as patches
21from obspy.io.sh.core import to_utcdatetime as toUTCDateTime
22from obspy.io.sh.core import from_utcdatetime as fromUTCDateTime
23
24DEG_TO_KM = 111.19
25
26provides = {"magnitude": "magnitude"}
27class magnitude(BaseCommand):
28    """
29    URI:http://www.seismic-handler.org/portal/wiki/ShMagnitude
30   
31    magnitude dump
32    magnitude dumpfile <magtype> <file>
33    magnitude determine <magtype> <trc> <from> <to>
34    magnitude determine <magtype> <trclist> <phaselist> <span>
35    magnitude mean <magtype> [<outvar>]
36    magnitude estimator <magtype> [<outvar-ml>] [<outvar-mlcorr>] [<picfile>]
37    magnitude clear <magtype> [<station>]
38    magnitude clear_all
39    magnitude set <magtype> <station> <value>
40    magnitude plot <magtype>
41
42    [Wiki Doc Text]
43
44    == Magnitude handling ==
45
46    '''command''': MAGNITUDE <subcmd> <magtype>[[br]]
47    or[[br]]
48    '''command''': MAGNITUDE DETERMINE <magtype> <trclist> <from> <to>[[br]]
49    or[[br]]
50    '''command''': MAGNITUDE DETERMINE <magtype> <trclist> <phaselist> <span>
51    or[[br]]
52    '''command''': MAGNITUDE SET <magtype> <station> <value>
53    or[[br]]
54    '''command''': MAGNITUDE PLOT <magtype> [<fname>]
55    or[[br]]
56    '''command''': MAGNITUDE DUMPFILE <magtype> <fname>
57
58    Magnitude management. Magnitudes are stored in a separate container, which is controlled via this command set.
59
60    === parameters ===
61
62     * <subcmd> ''parameter type: string''[[br]]
63       Magnitude subcommand: possible values are:[[br]]
64       DETERMINE: determine magnitude from waveform data[[br]]
65       SET: set magnitude values (mainly used for recovery operations)[[br]]
66       CLEAR: delete a subset of magnitudes[[br]]
67       CLEAR_ALL: delete all magnitudes[[br]]
68       MEAN: compute mean magnitude[[br]]
69       PLOT: create magnitude vs distance plot[[br]]
70       DUMP: dump all magnitudes on screen.
71       DUMPFILE: dump all magnitudes and distances of a magnitude type to file
72
73     * <magtype> ''parameter type: string''[[br]]
74       Magnitude type, possible values: 'ml', 'mb', 'mS', 'mBB'.
75
76     * <trclist> ''parameter type: trace list''[[br]]
77       List of traces to compute magnitudes.
78
79     * <from> ''parameter type: float''[[br]]
80       Start of (fixed) time window to find magnitude on waveforms.
81
82     * <to> ''parameter type: float''[[br]]
83       End of (fixed) time window to find magnitude on waveforms.
84
85     * <phaselist> ''parameter type: string''[[br]]
86       Comma separated list of phases. Find magnitude windows around phases. Uses manual (preferred) and theoretical phases. Takes the largest magnitude found in all windows.
87
88     * <span> ''parameter type: float''[[br]]
89       Width of time span to find magnitudes around phases. <span> is taken after manual phases or around (centered) theoretical phases.
90
91    === qualifiers ===
92
93    None
94
95    === example ===
96
97     `@MAGNITUDE DETERMINE ml ALL Pg,Sg 8.`::
98     Determine ml-magnitudes in 8s windows around Pg and Sg phases on all traces.
99
100     `magnitude mean ml &g1`::
101     Compute mean ml magnitude on all stations and put the result into variable g1.
102
103     `magnitude plot ml`::
104     Generate magnitude vs distance plot with all ml magnitudes and show figure on screen.
105
106     `magnitude clear ml`::
107     Delete all ml magnitudes.
108
109     `magnitude clear ml gr.gra1`::
110     Delete ml magnitude on station gr.gra1.
111    """
112    numberOfParameters = [1,2,3,4,5,6]
113    parameterQueries = [
114        {
115            "text": "subcmd",
116            "type": "str",
117            "question": False,
118        },
119        {
120            "text": "magtype",
121            "type": "str",
122            "question": False,
123        },
124        {
125            "text": "trace (list)",
126            "type": "str",
127            "question": False,
128        },
129        {
130            "text": "phaselist/from",
131            "type": "str",
132            "question": False,
133        },
134        {
135            "text": "span/to",
136            "type": "str",
137            "question": False,
138        },
139    ]
140    known_qualifiers = [
141    ]
142    redraw = False
143   
144    legal_subcmds = {
145        'dump'      : 1,
146        'dumpfile'  : 3,
147        'clear_all' : 1,
148        'plot'      : 2,
149        'determine' : 5,
150        'mean'      : 2,
151        'clear'     : 2,
152        'set'       : 4,
153        'estimator' : 3,
154    }
155   
156    def __init__(self, *args, **kwargs):
157        # unroll args & kwargs
158        BaseCommand.__init__(self, *args, **kwargs)
159
160    #@timeit
161    def run(self):
162        "Interface to magnitude management."
163
164        subcmd = self.parameters[0].lower()
165        if subcmd not in self.legal_subcmds.keys():
166            raise ShxError( "unknown magnitude subcmd '%s'" % subcmd )
167        if len(self.parameters) < self.legal_subcmds[subcmd]:
168            raise ShxError( "illegal number of parameters" )
169       
170        maglist = MagnitudeList()
171       
172        if subcmd == 'dump':
173            for station in maglist.getStations():
174                for mag in maglist.getMagnitudeTypes(station):
175                    print "%7s  %4s  %g" % (station,mag,
176                        maglist.getStationMagnitude(station,mag))
177        elif subcmd == 'mean':
178            magtype = self.getMagtype( self.parameters[1] )
179            mm = maglist.meanMagnitude( magtype )
180            if len(self.parameters) > 2:
181                magsymb = self.parameters[2][1:]
182                self.symbols.set( magsymb, mm )
183            else:
184                print mm
185        elif subcmd == 'clear_all':
186            maglist.clearAll()
187            set_infodisplay( "all magnitudes cleared" )
188        elif subcmd == 'clear':
189            magtype = self.getMagtype( self.parameters[1] )
190            if len(self.parameters) > 2:
191                station = self.parameters[2]
192                # convert trace number to station, if trace number given
193                try:
194                    trcno = int( station )
195                except:
196                    trcno == None
197                if trcno != None:
198                    try:
199                        trc = traces_from_list( "%d" % trcno )[0]
200                    except:
201                        return
202                    station = "%s.%s" % (trc.stats.network,trc.stats.station)
203                maglist.clearStationMagnitude( station, magtype )
204                set_infodisplay( "cleared %s on %s" % (magtype,station) )
205            else:
206                maglist.clearMagnitude( magtype )
207        elif subcmd == 'set':
208            magtype = self.getMagtype( self.parameters[1] )
209            station = self.parameters[2]
210            value = float( self.parameters[3] )
211            maglist.addStationMagnitude( station, magtype, value )
212        elif subcmd == 'plot':
213            magtype = self.getMagtype( self.parameters[1] )
214            if len(self.parameters) > 2:
215                picfile = self.parameters[2].lower()
216            else:
217                picfile = None
218            meanmag = maglist.meanMagnitude( magtype )
219            if meanmag == None:
220                title = "Magnitudes (%s) vs Distance (no mean)" % magtype
221            else:
222                title = "Magnitudes (%s) vs Distance (mean %4.2f)" % (magtype,
223                    meanmag)
224            fig = plt.figure( facecolor='white' )
225            #fig.suptitle( title )
226            ax = fig.add_subplot(111)
227            ax.set_xlabel( "distance (km)" )
228            ax.set_ylabel( "Magnitude" )
229            ax.set_title( title )
230            #ax.set_xlim( 0., 700. )
231            #ax.set_ylim( 0.1, 1000. )
232            dists, mags, labels = self.getMagnitudePlotData( magtype )
233            ax.scatter( dists, mags )
234            if min(mags) < 0.:
235                ax.set_ylim( 0., 1.1*max(mags) )
236            for lab,x,y in zip(labels,dists,mags):
237                plt.annotate( lab, xy=(x,y), fontsize=9 )
238            plt.draw()
239            showit = False
240            if picfile == None:
241                showit = True
242                picfile = self.tmpPlotfile()
243            plt.savefig( picfile, facecolor=fig.get_facecolor() )
244            if showit:
245                try:
246                    dspprog = Settings.config.misc.display_prog
247                except:
248                    dspprog = 'display'
249                os.system( "%s %s &" % (dspprog,picfile) )
250        elif subcmd == 'dumpfile':
251            magtype = self.getMagtype( self.parameters[1] )
252            dumpfile = self.parameters[2].lower()
253            if dumpfile.endswith('.png'):
254                dumpfile = os.path.splitext(dumpfile)[0] + '.txt'
255            fp = open( dumpfile, 'w' )
256            for x,y,lab in zip(*self.getMagnitudePlotData(magtype)):
257                fp.write( "%-8s %g %g\n" % (lab,x,y) )
258            fp.close()
259        elif subcmd == 'determine':
260            ap = AnalysisPar()
261            epilat = ap.getValue( 'epi_latitude' )
262            epilon = ap.getValue( 'epi_longitude' )
263            if epilat == None or epilon == None:
264                print "No epicenter -> no magnitude"
265                return
266            magtype = self.getMagtype( self.parameters[1] )
267            trclist = traces_from_list( self.parameters[2] )
268            reslist = {}
269            for trc in trclist:
270                # station name, skip if already in result list
271                sname = "%s.%s" % (trc.stats.network,trc.stats.station)
272                if sname in reslist:
273                    continue
274                # spans contains all time spans searched for magnitudes
275                spans = []
276                try:
277                    # if given explicitely
278                    stime = float( self.parameters[3] )
279                except:
280                    # given as phase list
281                    stime = None
282                etime = float( self.parameters[4] )
283                if stime == None:
284                    for phase in self.parameters[3].split(','):
285                        stime, etime = self.getPhaseWindow( trc, phase, etime )
286                        if stime != None:
287                            spans.append( (stime,etime) )
288                else:
289                    spans.append( (stime,etime) )
290                # determine maximum magnitude in all spans
291                degdist = getDistance( trc, epilat, epilon )
292                if degdist == None:
293                    print "No distance for '%s'" % sname
294                    continue
295                mmax = None
296                for stime,etime in spans:
297                    m = self.determineMagnitude( trc, magtype, stime, etime,
298                        degdist )
299                    if mmax == None or m > mmax:
300                        mmax = m
301                reslist[sname] = mmax
302            magcnt = 0
303            lastmag = None
304            laststation = None
305            for sname in reslist.keys():
306                if reslist[sname] != None:
307                    if magtype == 'mlz':
308                        lmagtype = 'ml'
309                    else:
310                        lmagtype = magtype
311                    maglist.addStationMagnitude( sname, lmagtype,
312                        reslist[sname] )
313                    magcnt += 1
314                    lastmag = reslist[sname]
315                    laststation = sname
316            if magcnt == 1:
317                set_infodisplay( "determined %s %4.2f on %s" % (magtype,lastmag,
318                    laststation) )
319            else:
320                set_infodisplay( "determined %d %s values" % (magcnt,magtype) )
321        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
328            magtype = self.getMagtype( self.parameters[1] )
329            if magtype != 'ml':
330                print "Estimator only for ml"
331                return
332            trclist = traces_from_list( self.parameters[2] )
333            reslist = []
334            stalist = []
335            for trc in trclist:
336                # station name, skip if already in result list
337                sname = "%s.%s" % (trc.stats.network,trc.stats.station)
338                if sname not in stalist:
339                    degdist = getDistance( trc, epilat, epilon )
340                    mag = maglist.getStationMagnitude(sname,magtype)
341                    if degdist != None and mag != None:
342                        reslist.append( (sname,degdist*DEG_TO_KM,mag) )
343                        stalist.append( sname )
344            me = MlEstimator()
345            magest, corrmag = me.estimateMl( reslist )
346            magest, corrmag = me.estimateMl( reslist, magnstart=magest )
347            if len(self.parameters) > 3:
348                magsymb = self.parameters[3][1:]
349                self.symbols.set( magsymb, magest )
350            else:
351                print "magn", magest
352            if len(self.parameters) > 4:
353                magsymb = self.parameters[4][1:]
354                self.symbols.set( magsymb, corrmag )
355            else:
356                print "corr magn", corrmag
357            if len(self.parameters) > 5:
358                picfile = self.parameters[5].lower()
359                me.plot( picfile )
360        else:
361            raise ShxError(
362                "program bug: magnitude subcmd '%s' not implemented" % subcmd )
363   
364    def getMagtype( self, mt ):
365        "Changes all uppercase names into correct magnitude notation."
366        mtl = mt.lower()
367        if mtl == 'ml':
368            return 'ml'
369        elif mtl == 'mlz':
370            return 'mlz'
371        elif mtl == 'mb':
372            return 'mb'
373        elif mtl == 'ms':
374            return 'mS'
375        elif mtl == 'mbb':
376            return 'mBB'
377        else:
378            return None
379   
380    def getPhaseWindow( self, trc, phase, wdwlen ):
381        pl = PhaseList()
382        manual = pl.translatePicktype( 'manual' )
383        auto = pl.translatePicktype( 'auto' )
384        theo = pl.translatePicktype( 'theo' )
385        sname = "%s.%s" % (trc.stats.network,trc.stats.station)
386        # first search manual picks
387        plist = pl.getPhaseList( sname, manual )
388        if len(plist) == 0:
389            plist = pl.getPhaseList( sname, auto )
390        for lphase in plist:
391            if lphase.name == phase:
392                reltime = lphase.picktime - trc.stats.starttime \
393                    - trc.get_info('t-origin')
394                if reltime == None:
395                    print "dbg: *** no reltime on", trc
396                    return
397                # manual phase -> window after pick
398                return (reltime,(reltime+wdwlen))
399        # then search theo picks
400        plist = pl.getPhaseList( sname, theo )
401        for lphase in plist:
402            if lphase.name == phase:
403                reltime = lphase.picktime - trc.stats.starttime \
404                    - trc.get_info('t-origin')
405                # theo phase -> window around pick
406                return ((reltime-wdwlen/2.),(reltime+wdwlen/2.))
407        return (None,None)
408   
409    def determineMagnitude( self, trc, magtype, stime, etime, degdist ):
410        if magtype == 'ml':
411            return self.determineMagnitudeMl( trc, stime, etime, degdist )
412        elif magtype == 'mlz':
413            return self.determineMagnitudeMlZ( trc, stime, etime, degdist )
414        else:
415            raise ShxError( "magnitude determination for '%s' not implmented"
416                % magtype )
417   
418    def determineMagnitudeMl( self, trc, stime, etime, degdist ):
419        # find horizontals
420        horizcomp = "NE12XY"
421        hlist = []
422        for t in traces_from_list('all'):
423            if t.stats.network != trc.stats.network:
424                continue
425            if t.stats.station != trc.stats.station:
426                continue
427            if t.stats.channel[-1].upper() not in horizcomp:
428                continue
429            # check simulation
430            if t.get_info('SIMUL') != "WOOD-ANDERSON":
431                continue
432            hlist.append( t )
433        if len(hlist) != 2:
434            print "No magnitude for %s: got %d traces" \
435                % (trc.stats.station,len(hlist))
436            return None
437        ampl = []
438        for t in hlist:
439            sl = trc.slice_relative( stime, etime )
440            #ampl.append( max( abs(max(sl.data)), abs(min(sl.data)) ) )
441            mean = sl.data.mean()
442            ampl.append( max(sl.data.max()-mean, mean-sl.data.min()) )
443        ampl = math.sqrt( ampl[0]*ampl[0] + ampl[1]*ampl[1] )
444        return magn_ml_analytic( ampl, degdist*DEG_TO_KM )
445   
446    def determineMagnitudeMlZ( self, trc, stime, etime, degdist ):
447        "Determine magnitude on Z component."
448        if trc.get_info('SIMUL') != "WOOD-ANDERSON":
449            print "No Wood-Anderson simulation"
450            return None
451        wdw = trc.slice_relative(stime,etime)
452        mean = wdw.data.mean()
453        wa_ampl_mean = max(wdw.data.max()-mean, mean-wdw.data.min())
454        h_dist = degdist * DEG_TO_KM
455        magnitude = math.log10(wa_ampl_mean) + math.log10(h_dist / 100.0) + \
456            0.00301 * (h_dist - 100.0) + 3.0
457        return magnitude
458   
459    def getMagnitudePlotData( self, magtype ):
460        maglist = MagnitudeList()
461        labels = []
462        dists = []
463        mags = []
464        ap = AnalysisPar()
465        epilat = ap.getValue( 'epi_latitude' )
466        epilon = ap.getValue( 'epi_longitude' )
467        for trc in traces_from_list('all'):
468            station = "%s.%s" % (trc.stats.network,trc.stats.station)
469            if station in labels:
470                continue
471            dist = getDistance( trc, epilat, epilon )
472            if dist == None:
473                continue
474            if magtype == 'ml':
475                dist *= DEG_TO_KM
476            m = maglist.getStationMagnitude( station, magtype )
477            if m == None:
478                continue
479            labels.append( station )
480            dists.append( dist )
481            mags.append( m )
482        return (dists,mags,labels)
483
484    def tmpPlotfile( self ):
485        return os.path.join( Settings.config.paths.scratch[0],
486            'tmpmagfig_%d.png' % os.getpid() )
487
488MAGN_ML_MIN_DIST = 1.0
489MAGN_ML_MAX_DIST = 1000.0
490
491def magn_ml_analytic( ampl, distance ):
492
493    """
494    /* Returns magnitude ml using analytical formula.  Amplitude and period must
495     * be read from a WOOD-ANDERSON simulated seismogram, horizontal component.
496     *
497     * parameters of routine
498     * float      ampl;        input; amplitude of WOOD-ANDERSON simulation
499     * float      distance;    input; epicentral distance in km
500     *                         returns value of ml
501     *
502     * Using analytical formula from IASPEI recommendations, see also
503     * http://www.seismic-handler.org/export/HEAD/SH_SHM/tags/2012a/doc/\
504     * Summary_WG-Recommendations_20110909.pdf
505         */
506    """
507    if distance < MAGN_ML_MIN_DIST or distance >= MAGN_ML_MAX_DIST:
508        return None
509    # analytical magnitude formula
510    ampl *= 1000 #/* nm */
511    magn = math.log10(ampl)
512    magn += 1.11 * math.log10(distance)
513    magn += 0.00189 * distance
514    magn -= 2.09
515    return magn
516
517class MlEstimator:
518
519    """Estimate ml magnitude of a preferably large set of station magnitudes.
520    Uses output file 'ml_statistics.at' of program ml_statistics.py.
521    """
522   
523    def __init__( self ):
524        self.initialized = False
525        self.plotdata = None
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
531        else:
532            self.initialized = False
533   
534    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
535
536    def parseFile( self, fname ):
537        self.statcorr = {}
538        for line in file(fname):
539            if line.startswith("#"):
540                continue
541            key = "constcorr:"
542            if line.startswith(key):
543                try:
544                    constcorr = float( line[len(key):] )
545                except:
546                    return False
547                if constcorr != 1:
548                    print "MlEstimator: need constant station corrections"
549                    return False
550                continue
551            key = "sensitivity range zeroval:"
552            if line.startswith(key):
553                try:
554                    self.rangestart = float( line[len(key):] )
555                except:
556                    print "Cannot read '%s'" % key
557                    return False
558                continue
559            key = "sensitivity range inc:"
560            if line.startswith(key):
561                try:
562                    self.rangeinc = float( line[len(key):] )
563                except:
564                    print "Cannot read '%s'" % key
565                    return False
566                continue
567            key = "2nd stage zeroval:"
568            if line.startswith(key):
569                try:
570                    self.secondordercorr_zero = float( line[len(key):] )
571                except:
572                    print "Cannot read '%s'" % key
573                    return False
574                continue
575            key = "2nd stage inc:"
576            if line.startswith(key):
577                try:
578                    self.secondordercorr_inc = float( line[len(key):] )
579                except:
580                    print "Cannot read '%s'" % key
581                    return False
582                continue
583            if ':' in line:
584                try:
585                    sstr, cstr = line.split(':')
586                    self.statcorr[sstr.strip()] = float( cstr.strip() )
587                except:
588                    print "parse error in line:", line
589                    return False
590        return True
591   
592    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
593
594    def estimateMl( self, values, magnstart=None ):
595        """Needs list of triples (station,distance,ml), station like 'GR.BFO',
596        distance in deg, ml is uncorrected station magnitude.
597        magnstart is a first approximation of the magnitude. If not specified
598        derived from the closest stations.
599        """
600        if magnstart == None:
601            magnstart = self.magnApprox( values )
602            if magnstart == None:
603                return (None,None)
604        mindist = 30.
605        maxdist = self.rangestart + (magnstart-1.)*self.rangeinc
606        dv = []
607        mv = []
608        for sta,dist,magn in values:
609            if dist < mindist or dist > maxdist:
610                continue
611            dv.append( dist )
612            if sta in self.statcorr.keys():
613                mv.append( magn-self.statcorr[sta] )
614            else:
615                mv.append( magn )
616        dmperkm, zeromag = np.polyfit( np.array(dv), np.array(mv), 1 )
617        corrmag = zeromag-self.secondOrderCorrection(zeromag,dmperkm)
618        self.plotdata = (values,mindist,maxdist,zeromag,dmperkm,corrmag)
619        return (zeromag,corrmag)
620           
621    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
622
623    def plot( self, fname ):
624        if self.plotdata == None:
625            return
626        values,mindist,maxdist,zeromag,dmperkm,corrmag = self.plotdata
627        fig = plt.figure( facecolor='white' )
628        ax = fig.add_subplot(111)
629        ax.set_xlabel( "distance (km)" )
630        ax.set_ylabel( "Magnitude" )
631        ax.set_title( "ml: %4.2f  corr ml: %4.2f" % (zeromag,corrmag) )
632        ax.set_xlim( 0., 500. )
633        ax.set_ylim( 0.4, 4.5 )
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 )
642        maxd = max( dists )
643        ax.add_patch(
644            patches.Rectangle(
645                (0.,0.),mindist,5.,facecolor='yellow',
646                edgecolor='white',alpha=0.5
647            )
648        )
649        ax.add_patch(
650            patches.Rectangle(
651                (maxdist,0.),500.-maxdist,5.,facecolor='yellow',
652                edgecolor='white',alpha=0.5
653            )
654        )
655        ax.plot( [0.,maxd], [zeromag,(zeromag+dmperkm*maxd)], c='r' )
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):
661            plt.annotate( lab, xy=(x,y), fontsize=9 )
662        plt.draw()
663        plt.savefig( fname )
664        plt.close( fig )
665       
666
667    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
668
669    def magnApprox( self, val ):
670        "Estimate magnitude from closest 10% of stations."
671        if len(val) < 20:
672            return None
673        sta, dist, mag = zip(*val)
674        limdist = sorted(dist)[len(dist)/10]
675        mags = []
676        for sta,dist,mag in val:
677            if dist <= limdist:
678                mags.append( mag )
679        if len(mags) > 4:
680            mags.remove( min(mags) )
681            mags.remove( max(mags) )
682        return (sum(mags)/float(len(mags)))
683
684    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
685
686    def secondOrderCorrection( self, magn, dmperkm ):
687        "Returns 2nd order correction for magnitude."
688        if self.secondordercorr_zero == None:
689            return 0.
690        else:
691            return (self.secondordercorr_zero \
692                + self.secondordercorr_inc*dmperkm)
693
Note: See TracBrowser for help on using the repository browser.