From a6d62489ee0f909a8c1dc5ec6366d11e37b02e55 Mon Sep 17 00:00:00 2001 From: Michael Hearne Date: Mon, 8 Apr 2013 14:54:23 -0600 Subject: [PATCH] Added support for focal mechanism angles at command line, and started working on code to support retrieval of MT info from USGS ComCat --- .gitignore | 4 ++ strec.py | 43 ++++++++++++- strec/gmpe.py | 4 +- strec/mtreader.py | 161 ++++++++++++++++++++++++++++++++++++++++++++-- strec_convert.py | 12 +++- strec_init.py | 7 ++ 6 files changed, 222 insertions(+), 9 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b502a9b --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*.py[cod] +*.grd +*.db +run* diff --git a/strec.py b/strec.py index 71086a0..24b1b68 100755 --- a/strec.py +++ b/strec.py @@ -6,9 +6,31 @@ import ConfigParser import getpass import datetime +import math from strec.gmpe import GMPESelector import strec.utils +from strec import cmt + +def getPlungeValues(strike,dip,rake,mag): + mom = 10**((mag*1.5)+16.1) + d2r = math.pi/180.0 + + mrr=mom*math.sin(2*dip*d2r)*math.sin(rake*d2r) + mtt=-mom*((math.sin(dip*d2r)*math.cos(rake*d2r)*math.sin(2*strike*d2r))+(math.sin(2*dip*d2r)*math.sin(rake*d2r)*(math.sin(strike*d2r)*math.sin(strike*d2r)))) + mpp=mom*((math.sin(dip*d2r)*math.cos(rake*d2r)*math.sin(2*strike*d2r))-(math.sin(2*dip*d2r)*math.sin(rake*d2r)*(math.cos(strike*d2r)*math.cos(strike*d2r)))) + mrt=-mom*((math.cos(dip*d2r)*math.cos(rake*d2r)*math.cos(strike*d2r))+(math.cos(2*dip*d2r)*math.sin(rake*d2r)*math.sin(strike*d2r))) + mrp=mom*((math.cos(dip*d2r)*math.cos(rake*d2r)*math.sin(strike*d2r))-(math.cos(2*dip*d2r)*math.sin(rake*d2r)*math.cos(strike*d2r))) + mtp=-mom*((math.sin(dip*d2r)*math.cos(rake*d2r)*math.cos(2*strike*d2r))+(0.5*math.sin(2*dip*d2r)*math.sin(rake*d2r)*math.sin(2*strike*d2r))) + + plungetuple = cmt.compToAxes(mrr,mtt,mpp,mrt,mrp,mtp) + plungevals = {} + plungevals['T'] = plungetuple[0].copy() + plungevals['N'] = plungetuple[1].copy() + plungevals['P'] = plungetuple[2].copy() + plungevals['NP1'] = plungetuple[3].copy() + plungevals['NP2'] = plungetuple[4].copy() + return plungevals if __name__ == '__main__': @@ -26,6 +48,9 @@ parser.add_option("-d", "--datafile",dest="datafile", metavar="DATAFILE", help="Specify the database (.db) file containing moment tensor solutions.") + parser.add_option("-a", "--angles",dest="angles", + metavar="ANGLES", + help='Specify the focal mechanism by providing "strike dip rake"') parser.add_option("-c", "--csv-out", action="store_true", dest="outputCSV", default=False, help="print output as csv") @@ -80,14 +105,30 @@ etime = None forceComposite = True + if options.angles is not None: + parts = options.angles.split() + try: + strike = float(parts[0]) + dip = float(parts[1]) + rake = float(parts[2]) + plungevals = getPlungeValues(strike,dip,rake,magnitude) + except: + print 'Could not parse Strike, Dip and Rake from "%s"' % options.angles + sys.exit(1) + else: + plungevals = None + if options.forceComposite: forceComposite = True + plungevals = None homedir = os.path.dirname(os.path.abspath(__file__)) #where is this script? zoneconfigfile = os.path.join(homedir,'strec.ini') gs = GMPESelector(zoneconfigfile,datafile,homedir,datafolder) - strecresults = gs.selectGMPE(lat,lon,depth,magnitude,date=etime,forceComposite=forceComposite) + strecresults = gs.selectGMPE(lat,lon,depth,magnitude,date=etime, + forceComposite=forceComposite, + plungevals=plungevals) if options.outputCSV: strecresults.renderCSV(sys.stdout) diff --git a/strec/gmpe.py b/strec/gmpe.py index 6b5410c..673bee8 100755 --- a/strec/gmpe.py +++ b/strec/gmpe.py @@ -267,7 +267,9 @@ def selectGMPE(self,lat,lon,depth,magnitude,date=None,forceComposite=False,plung plungevals,fmwarning,mtsource = self.getFocMechAxes(lat,lon,depth,magnitude) else: plungevals,fmwarning,mtsource = self.getFocMechAxes(lat,lon,depth,magnitude,date=date) - + else: + mtsource = 'input' + fmwarning = '' if plungevals is None: plungevals = {'T':{},'N':{},'P':{},'NP1':{},'NP2':{}} plungevals['T']['plunge'] = numpy.nan diff --git a/strec/mtreader.py b/strec/mtreader.py index 6ce4e1c..8cbc916 100755 --- a/strec/mtreader.py +++ b/strec/mtreader.py @@ -1,15 +1,24 @@ #!/usr/bin/env python +#stdlib imports import re import datetime +import pytz import math from xml.dom.minidom import parse -from cmt import compToAxes import os.path import sqlite3 +import urllib,urllib2 +import json +import sys +import copy + +#local imports +from cmt import compToAxes + -class CSVReaderError(Exception): - """Used to indicate errors in the CSV reader class.""" +class ReaderError(Exception): + """Used to indicate errors in the reader classes.""" def appendDataFile(infile,dbfile,filetype,dtype,hasHeader=False): if filetype == 'ndk': @@ -84,6 +93,135 @@ def __init__(self,mtfile,type=None): def generateRecords(self,startDate=None,enddate=None,hasHeader=False): pass +class ComCatReader(MTReader): + URLBASE = 'http://comcat.cr.usgs.gov/earthquakes/feed/search.php?%s' + def __init__(self): + pass + + def getEvents(self,pdict): + params = urllib.urlencode(pdict) + searchurl = self.URLBASE % params + datadict = None + try: + fh = urllib2.urlopen(searchurl) + data = fh.read() + data2 = data[len(pdict['callback'])+1:-2] + datadict = json.loads(data2) + fh.close() + except Exception,errorobj: + raise errorobj + return datadict['features'] + + def generateRecords(self,startdate=None,enddate=None,hasHeader=False): + utc = pytz.UTC + if startdate is None: + startdate = datetime.datetime(1971,1,1) + startdate = utc.localize(startdate) + if enddate is None: + enddate = datetime.datetime(1971,1,1) + enddate = utc.localize(enddate) + if (enddate-startdate).days < 365: + mintime = int(startdate.strftime('%s'))*1000 + maxtime = int(enddate.strftime('%s'))*1000 + params = pdict = {'productType[]':'moment-tensor,', + 'minEventTime':mintime, + 'maxEventTime':maxtime, + 'callback':'comsearch'} + events = self.getEvents(params) + else: + ndays = (enddate-startdate).days + nyears = int(math.ceil(ndays/365.0)) + events = [] + mindate = startdate + maxdate = mindate + datetime.timedelta(days=365) + for i in range(0,nyears): + mintime = int(mindate.strftime('%s'))*1000 + maxtime = int(maxdate.strftime('%s'))*1000 + params = pdict = {'productType[]':'moment-tensor,', + 'minEventTime':mintime, + 'maxEventTime':maxtime, + 'callback':'comsearch'} + events = events + self.getEvents(params) + mindate = maxdate + maxdate = mindate + datetime.timedelta(days=365) + for event in events: + etimesecs = int(event['properties']['time'])/1000 + etime = datetime.datetime.utcfromtimestamp(etimesecs) + if etime < startdate or etime > enddate: + continue + edicts = self.getMomentDicts(event) + for edict in edicts: + yield edict + + def getMomentDicts(self,event): + eventdicts = [] + datadict = {} + url = event['properties']['url'] + '.json' + try: + fh = urllib2.urlopen(url) + data = fh.read() + datadict = json.loads(data) + fh.close() + except Exception,errorobj: + raise errorobj + eqtime = datadict['summary']['time'] + eqid = str(datetime.datetime.utcfromtimestamp(int(eqtime)/1000)) + lat = float(datadict['summary']['latitude']) + lon = float(datadict['summary']['longitude']) + depth = float(datadict['summary']['depth']) + mag = float(datadict['summary']['magnitude']) + for mt in datadict['products']['moment-tensor']: + if mt['status'] == 'DELETE': + continue + bsource = mt['properties']['beachball-source'] + bcode = mt['code'] + #we have gcmt data going back further, directly from GCMT website + if bsource == 'LD' and bcode.find('GCMT') > -1: + continue + + mttype = mt['properties']['beachball-type'] + try: + mrr = float(mt['properties']['tensor-mrr']) + mtt = float(mt['properties']['tensor-mtt']) + mpp = float(mt['properties']['tensor-mpp']) + mrt = float(mt['properties']['tensor-mrt']) + mrp = float(mt['properties']['tensor-mrp']) + mtp = float(mt['properties']['tensor-mtp']) + except: + continue + np1strike = float(mt['properties']['nodal-plane-1-strike']) + np2strike = float(mt['properties']['nodal-plane-2-strike']) + np1dip = float(mt['properties']['nodal-plane-1-dip']) + np2dip = float(mt['properties']['nodal-plane-2-dip']) + try: + np1rake = float(mt['properties']['nodal-plane-1-rake']) + except: + np1rake = float(mt['properties']['nodal-plane-1-slip']) + try: + np2rake = float(mt['properties']['nodal-plane-2-rake']) + except: + np2rake = float(mt['properties']['nodal-plane-2-slip']) + try: + T,N,P,NP1,NP2 = compToAxes(mrr,mtt,mpp,mrt,mrp,mtp) + except Exception,exc: + pass + tplunge = T['plunge'] + tazimuth = T['azimuth'] + nplunge = N['plunge'] + nazimuth = N['azimuth'] + pplunge = P['plunge'] + pazimuth = P['azimuth'] + eventdicts.append({'id':eqid,'type':mttype, + 'lat':lat,'lon':lon, + 'depth':depth,'mag':mag, + 'mrr':mrr,'mtt':mtt, + 'mpp':mpp,'mrt':mrt, + 'mrp':mrp,'mtp':mtp, + 'tplunge':tplunge,'tazimuth':tazimuth, + 'nplunge':nplunge,'nazimuth':nazimuth, + 'pplunge':pplunge,'pazimuth':pazimuth}) + return eventdicts + class CSVReader(MTReader): idfmtsec = '%Y%m%d%H%M%S' idfmtmin = '%Y%m%d%H%M' @@ -102,7 +240,7 @@ def generateRecords(self,startdate=None,enddate=None,hasHeader=False): else: yield self.readline(line) except Exception,msg: - raise CSVReaderError,"Error reading line %s" % line + raise ReaderError,'Error reading data from line:\n"%s".\nIs this a header row?' % line.strip() self.fh.close() @@ -409,10 +547,23 @@ def generateRecords(self,startdate=None,enddate=None): self.dom.unlink() if __name__ == '__main__': + #get a months worth of moment tensor records from ComCat + comreader = ComCatReader() + d1 = datetime.datetime(1976,1,1,0,0) + #d2 = datetime.datetime.now() + #d2 = datetime.datetime(2002,10,1) + d2 = datetime.datetime.now() + nc = 0 + for record in comreader.generateRecords(startdate=d1,enddate=d2): + print '%s,%.1f,%s' % (record['id'],record['mag'],record['type']) + sys.stdout.flush() + nc += 1 + print 'There are %i moment tensor records in ComCat' % nc + sys.exit(0) xreader = QuakeMLReader('xmlsample.xml',type='GCMT') for record in xreader.generateRecords(): print '%s,%.1f,%s' % (record['id'],record['mag'],record['type']) - + creader = CSVReader('csvtest.csv','User') for record in creader.generateRecords(): print record['id'],record['mag'],record['type'] diff --git a/strec_convert.py b/strec_convert.py index b039ab3..be15e09 100755 --- a/strec_convert.py +++ b/strec_convert.py @@ -82,6 +82,14 @@ if os.path.isfile(outfile): print '%s already exists - appending new data.' % outfile - appendDataFile(infile,outfile,filetype,options.fmtype,hasHeader=options.hasheader) + try: + appendDataFile(infile,outfile,filetype,options.fmtype,hasHeader=options.hasheader) + except Exception,msg: + print 'Error reading input file %s.\n%s' % (infile,msg) + sys.exit(1) else: - createDataFile(infile,outfile,filetype,options.fmtype,hasHeader=options.hasheader) + try: + createDataFile(infile,outfile,filetype,options.fmtype,hasHeader=options.hasheader) + except Exception,msg: + print 'Error reading input file %s.\n"%s"' % (infile,msg) + sys.exit(1) diff --git a/strec_init.py b/strec_init.py index 35b9e8e..a90fef6 100755 --- a/strec_init.py +++ b/strec_init.py @@ -128,6 +128,9 @@ def fetchSlabs(datafolder): parser.add_option("-g", "--gcmt", action="store_true", dest="getGCMT", default=False, help="Download all GCMT moment tensor data") + parser.add_option("-c", "--comcat", + action="store_true", dest="getComCat", default=False, + help="Download all USGS ComCat moment tensor data (sans GCMT)") parser.add_option("-n", "--noslab", action="store_true", dest="noSlab", default=False, help="Do NOT download slab data") @@ -196,6 +199,9 @@ def fetchSlabs(datafolder): createDataFile(histfile,outfile,'ndk','gcmt',hasHeader=False) print 'Finished converting historical GCMT data.' os.remove(histfile) + if options.getComCat: + + if options.update: mostrecent = getMostRecent(outfile) ryear = mostrecent.year @@ -226,6 +232,7 @@ def fetchSlabs(datafolder): mostrecent = getMostRecent(outfile) print 'GCMT database contains events through %s.' % str(mostrecent) sys.exit(0) +