 Timestamp:
 12/19/16 14:38:41 (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

SHX/trunk/SeismicHandler/commands/magnitude.py
r1204 r1206 7 7 import os 8 8 import math 9 import numpy as np 9 10 from SeismicHandler.basics.command import BaseCommand 10 11 from SeismicHandler.basics.analysispar import PhaseList, \ … … 32 33 magnitude determine <magtype> <trclist> <phaselist> <span> 33 34 magnitude mean <magtype> [<outvar>] 35 magnitude estimator <magtype> [<outvarml>] [<outvarmlcorr>] [<picfile>] 34 36 magnitude clear <magtype> [<station>] 35 37 magnitude clear_all … … 315 317 else: 316 318 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 ) 317 350 else: 318 351 raise ShxError( … … 472 505 return magn 473 506 507 class 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 + (magnstart1.)*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( magnself.statcorr[sta] ) 603 else: 604 mv.append( magn ) 605 dmperkm, zeromag = np.polyfit( np.array(dv), np.array(mv), 1 ) 606 corrmag = zeromagself.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.