Ground Motion Prediction¶
Import pygmm module¶
In [1]:
import pygmm
Import ObsPy module¶
In [2]:
from obspy import read
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.geodetics.base import gps2dist_azimuth
import obspy as ob
print("# obspy version = ",ob.__version__)
Import libcomcat module¶
In [3]:
# Local imports
from libcomcat.dataframes import get_phase_dataframe, get_magnitude_data_frame
from libcomcat.search import get_event_by_id
from libcomcat.dataframes import (get_detail_data_frame, get_dyfi_data_frame,
get_history_data_frame, get_magnitude_data_frame,
get_pager_data_frame, get_phase_dataframe,
get_summary_data_frame)
Import SciPy, NumPy, matplotlib module¶
In [4]:
import numpy as np
import scipy as sp
import matplotlib as mpl
import pandas as pd
print("# numpy version = ",np.__version__)
print("# scipy version = ",sp.__version__)
print("# matplotlib version = ",mpl.__version__)
print("# pandas version = ",pd.__version__)
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.dates as mdates
from scipy import signal
from scipy import ndimage
from scipy.fftpack import fft, ifft
from scipy.linalg import norm
from scipy.stats import pearsonr
from scipy.stats import spearmanr
from scipy.stats import kendalltau
import sys
import os
Font size for plotting¶
In [5]:
# font size
SMALL_SIZE = 16
MEDIUM_SIZE = 18
BIGGER_SIZE = 20
#SMALL_SIZE = 32
#MEDIUM_SIZE = 32
#BIGGER_SIZE = 36
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
Input parameters¶
In [6]:
# event id
eventid = "nc72948801" # event_id can be obtained from the spreadsheet
# target location
stla = 37.87352 # latitude of the target area # this example use BK.BRK station
stlo = -122.26099 # longitude of the target area # this example use BK.BRK station
# vs30
v_s30 = 500 # m/s find this vs30 from grd file
# https://github.com/usgs/earthquake-global_vs30/tree/master/California/California_Vs30_7p5c.grd
# region.
region = "california" # need for BSSA14
#mechanism="U" # this part need to be updated
Search for earthquake¶
This example use USGS event catalog with eventid
In [7]:
clientEQ = Client("USGS")
catalog = clientEQ.get_events(eventid=eventid)
print(catalog)
_plot = catalog.plot()
Extract event information¶
In [8]:
# event info. origin time, location, magnitude
event = catalog[0]
origin = event.origins[0]
#origin_time = origin.time
evla = origin.latitude
evlo = origin.longitude
evdp_km = origin.depth / 1000
evmag = event.magnitudes[0].mag
#evyearOUT = origin_time.year
#evjdayOUT = origin_time.julday
#evhourOUT = origin_time.hour
#evminOUT = origin_time.minute
#evsecOUT = origin_time.second
#evid = event.origins[0]['extra']['dataid']['value']
#event_region = event.event_descriptions[0]['text']
#print("# evid = ", evid)
Extract focal mechanism information¶
In [9]:
eq_detail = get_event_by_id(eventid, includesuperseded=True)
In [10]:
print("Included products: %s " % eq_detail.products)
product = eq_detail.getProducts('focal-mechanism')[0]
dip = product['nodal-plane-1-dip']
rake = product['nodal-plane-1-rake']
strike = product['nodal-plane-1-strike']
print("# strike = ", strike)
print("# dip = ", dip)
print("# rake = ", rake)
Create CSV file for "mechanism" estimate¶
output file is focal.csv
In [11]:
### use FMC.py to get fault-type
#AR Focal mechanism one plane (psmeca compatible):[longitude, latitude, depth, strike, dip, rake, magnitude (Mw),X plot, Y plot (for GMT), ID, TYPE]
focal_csv = "focal.csv"
#focal_df = pd.DataFrame({ 'lon' : evlo})
focal_df = pd.DataFrame({ 'lon' : evlo,
'lat' : evla,
'dep' : evdp_km,
'strike' : strike,
'dip' : dip,
'rake' : rake,
'mag' : evmag,
'evid' : eventid }, index=[0])
#focal_data = pd.read_csv(hash_catalog_hypoDD_fi ,
# sep=" ",names=["evid", "evyear", "evmonth", "evday", "evhour", "evmin", "evsec", "evtype", "mag", "magtype",
# "lat", "lon", "dep", "loc_quality", "RMS", "err_h", "err_d", "err_t", "num_pick_loc", "num_Ppick", "num_Spick",
# "strike", "dip", "rake", "fault_plane_uncert", "aux_plane_uncert", "num_p_pola", "wei_per_misfit_firstM", "focal_mech_qual", "prob_mech_sol", "sta_dist_ratio"], header=None)
#print(hash_catalog_hypoDD_data)
focal_df.to_csv(focal_csv, columns=['lon','lat','dep','strike','dip','rake','mag','lon','lat','evid'], header=False, index=False, sep=' ')
Compute "Rupture_type" with FMC.py¶
https://github.com/ElsevierSoftwareX/SOFTX_2018_227
output are focal.csv.png and focal.csv_fault_type.out
In [12]:
%%script tcsh
./SOFTX_2018_227-master/FMC.py < ./focal.csv -i AR -o ALL -p focal.csv.png >! focal.csv_fault_type.out
Read focal.csv_fault_type.out to get "Rupture_type" and convet it to "mechanism" ¶¶
In [13]:
#20920 N
#14252 N-SS
#13577 R
#12723 R-SS
#39297 SS
#30421 SS-N
#27145 SS-R
fault_type_fi = "focal.csv_fault_type.out"
fault_type_df = pd.read_csv(fault_type_fi,
sep=" ")
#print(fpfit_catalog_hypoDD_fault_type )
#Abbreviation Name
#u Unspecified
#ss Strike-slip
#ns Normal slip
#rs Reverse slip
fault_type_df ['fault_type'] = -9
fault_type_df ['mechanism'] = "U"
fault_type_df.loc[ (fault_type_df['Rupture_type'] == "N") | (fault_type_df['Rupture_type'] == "N-SS") , 'fault_type'] = 0
fault_type_df.loc[ (fault_type_df['Rupture_type'] == "N") | (fault_type_df['Rupture_type'] == "N-SS") , 'mechanism'] = "NS"
fault_type_df.loc[ (fault_type_df['Rupture_type'] == "R") | (fault_type_df['Rupture_type'] == "R-SS") , 'fault_type'] = 1
fault_type_df.loc[ (fault_type_df['Rupture_type'] == "R") | (fault_type_df['Rupture_type'] == "R-SS") , 'mechanism'] = "RS"
fault_type_df.loc[ (fault_type_df['Rupture_type'] == "SS") | (fault_type_df['Rupture_type'] == "SS-N") | (fault_type_df['Rupture_type'] == "SS-R") , 'fault_type'] = 2
fault_type_df.loc[ (fault_type_df['Rupture_type'] == "SS") | (fault_type_df['Rupture_type'] == "SS-N") | (fault_type_df['Rupture_type'] == "SS-R") , 'mechanism'] = "SS"
# to see if we have -9; if so, this would be a problem
#fault_type_df[fault_type_df['fault_type'] == -9]
Check Rupture_type and mechanism¶
In [14]:
print("# Rupture_type = ", fault_type_df['Rupture_type'])
In [15]:
mechanism = fault_type_df['mechanism'][0]
print("# mechanism = ", mechanism)
Computer JB distance¶
This example uses hypocenteral distance between the epicenter and terget area
In [16]:
#baz = gps2dist_azimuth(source_latitude, source_longitude, station_latitude, station_longitude)
baz = gps2dist_azimuth(evla, evlo, stla, stlo)
#print('Epicentral distance [m]: ', baz[0])
#print('Theoretical azimuth [deg]: ', baz[1])
#print('Theoretical backazimuth [deg]: ', baz[2])
dist_jb = baz[0]/1000 # convert to m -> km
print("# dist_jb = ", dist_jb)
Check input parameter for ground motion prediction¶
In [17]:
print("# evmag = ", evmag)
print("# dist_jb = ", dist_jb)
print("# v_s30 = ", v_s30)
print("# region = ", region)
print("# mechanism = ", mechanism)
Compute ground motion with BSSA14 model¶
In [18]:
scenario = pygmm.Scenario(mag=evmag, dist_jb=dist_jb, v_s30=v_s30, region=region, mechanism=mechanism)
m = pygmm.BooreStewartSeyhanAtkinson2014(scenario)
print("# PGA (g) = ", m.pga)
print("# PGV (cm/s) = ", m.pgv)
In [ ]:
In [ ]:
In [ ]:
In [ ]: