Skip to content
This repository has been archived by the owner on Jul 14, 2022. It is now read-only.

Commit

Permalink
Merge pull request #49 from emthompson-usgs/update
Browse files Browse the repository at this point in the history
New module for doing event-specific point source adjustments
  • Loading branch information
cbworden authored Jan 23, 2019
2 parents 45fe3a9 + 169f2cc commit 19ebf28
Show file tree
Hide file tree
Showing 7 changed files with 345 additions and 179 deletions.
3 changes: 2 additions & 1 deletion bin/run_ps2ff_single_event
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ The configuration parameters are:
*Seismological Research Letters*, 81(6), 941--950.
- Sea10_slab: Slab coefficients from the paper in previous bullet.
- *mech* -- The rupture mechanism, only used by some scaling relationships:
- **mech** -- The rupture mechanism, only used by some scaling relationships:
- A -- all/unknown mechanisms,
- SS -- strike-slip,
Expand Down Expand Up @@ -121,6 +121,7 @@ configuration file for this program is given in
help='Print informational messages.')
return parser


def main(args):

start_datetime = datetime.datetime.now().isoformat()
Expand Down
39 changes: 39 additions & 0 deletions code.json
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,45 @@
"measurementType": {},
"agency": "U.S. Department of Interior",
"releases": [
{
"name": "ps2ff",
"description": "Approximated rupture distances from point source",
"version": "1.3",
"status": "Production",
"organization": "U.S. Geological Survey",
"vcs": "git",
"repositoryURL": "https://github.com/usgs/ps2ff",
"homepageURL": "https://github.com/usgs/ps2ff",
"downloadURL": "https://github.com/usgs/ps2ff/releases/tag/v1.3",
"disclaimerURL": "https://github.com/usgs/ps2ff/blob/master/DISCLAIMER.md",
"permissions": {
"licenses": [
{
"name": "Public Domain",
"URL": "https://github.com/usgs/ps2ff/blob/master/LICENSE.md"
}
],
"usageType": "openSource",
"exemptionText": null
},
"laborHours": 1,
"languages": [ "Python" ],
"tags": [
"earthquake",
"shakemap",
"point source",
"finite fault"
"seismic-hazard",
"earthquake-hazard"
],
"contact": {
"name": "Eric Thompson",
"email": "[email protected]"
},
"date": {
"metadataLastUpdated": "2018-12-12"
}
},
{
"name": "ps2ff",
"description": "Approximated rupture distances from point source",
Expand Down
21 changes: 11 additions & 10 deletions ps2ff/integration_loops.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,9 +49,9 @@ def mag_dist_loop(conf, iP=None, filename=None, M=None, Repi=None):
ratio[i, j] = dist_avg / Repi[i]
variance[i, j] = dist_var
print(" Proc %d j=%d of %d %s" % (iP, j+1, nmag,
datetime.datetime.now().isoformat()))
datetime.datetime.now().isoformat()))
print("Proc %d done with %d of %d distances %s" % (iP, i+1, nepi,
datetime.datetime.now().isoformat()))
datetime.datetime.now().isoformat()))

for fname, farr in (('Ratios', ratio), ('Var', variance)):
f = open('%s%s_%02d.csv' % (filename, fname, iP), 'w')
Expand Down Expand Up @@ -110,7 +110,7 @@ def single_event_loop(conf, iP=0, rjb_filename='junk',
Rjb_ratio[i, j] = rjb_avg / Repi[i]
Rjb_variance[i, j] = rjb_var
print(" Proc %d j=%d of %d %s" % (iP, j+1,
conf['ntheta'], datetime.datetime.now().isoformat()))
conf['ntheta'], datetime.datetime.now().isoformat()))
print("Proc %d done with %d of %d distances %s" %
(iP, i+1, nepi, datetime.datetime.now().isoformat()))
else:
Expand All @@ -122,9 +122,9 @@ def single_event_loop(conf, iP=0, rjb_filename='junk',
Rjb_ratio[i, 0] = rjb_avg / Repi[i]
Rjb_variance[i, 0] = rjb_var
print(" Proc %d j=%d of %d %s" % (iP, 1, 1,
datetime.datetime.now().isoformat()))
datetime.datetime.now().isoformat()))
print("Proc %d done with %d of %d distances %s" % (iP, i+1, nepi,
datetime.datetime.now().isoformat()))
datetime.datetime.now().isoformat()))
fr_rrup = open('%sRatios_%02d.csv' % (rrup_filename, iP), 'w')
fv_rrup = open('%sVar_%02d.csv' % (rrup_filename, iP), 'w')
fr_rjb = open('%sRatios_%02d.csv' % (rjb_filename, iP), 'w')
Expand Down Expand Up @@ -167,7 +167,7 @@ def rjb_inner_loop(M, Repi, conf):
"""
This function evaluates the Rjb mean and var integral
for a single M/R pair, looping over:
- dip
- dx, dy (location of hypocenter on fault)
- theta (angle to fault)
Expand Down Expand Up @@ -504,7 +504,8 @@ def rrup_inner_loop(M, Repi, conf):
dx = x[1] - x[0]
# Calclate range of Ztor
ZtorMax = conf['max_seis_depth'] - width[w]*np.sin(dip[k])
if np.allclose(ZtorMax, 0) or conf['nz'] == 1: # Should 0 be min_seis_depth
# Should 0 be min_seis_depth
if np.allclose(ZtorMax, 0) or conf['nz'] == 1:
nz = 1
dz = 0
Ztor = np.array([0]) # should be min_seis_depth ???
Expand Down Expand Up @@ -597,8 +598,8 @@ def rrup_inner_loop(M, Repi, conf):
Ztor[z]*np.cos(dip[k])
r3 = Rx > (tmp + width[w]/np.cos(dip[k]))
Rrupp[r3] = np.sqrt((Rx[r3] -
width[w]*np.cos(dip[k]))**2 + (Ztor[z] +
width[w]*np.sin(dip[k]))**2)
width[w]*np.cos(dip[k]))**2 + (Ztor[z] +
width[w]*np.sin(dip[k]))**2)

Rrup = np.sqrt(Rrupp**2 + Ry**2)
Rrup2 = Rrup * Rrup
Expand Down Expand Up @@ -662,7 +663,7 @@ def single_event_inner_loop(conf, Repi, theta=0, ntheta=73):
Repi (float): Epicentral distance (km).
theta (float): Source-to-site angle (radians).
ntheta (int): Number of integration steps for theta; used if `bytheta`
is True.
is False.
Returns:
tuple: Rrup variance, Rrup mean, Rjb variance, Rjb mean
Expand Down
126 changes: 126 additions & 0 deletions ps2ff/run.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,126 @@

import numpy as np
import time as time

from ps2ff.integration_loops import single_event_inner_loop
from ps2ff.constants import MagScaling, Mechanism


# Use max_workers threading in model?

def single_event_adjustment(
magnitude,
hyp_depth,
ar=1.7,
mechanism=Mechanism.A,
mag_scaling=MagScaling.WC94,
n_repi=13,
min_repi=0.1,
max_repi=1000,
nxny=3,
n_theta=9,
n_dip=3,
min_dip=0,
max_dip=np.pi/2.0,
n_eps=3,
trunc=2):
"""
Method for getting event-specific point source distance
adjustment factors. This does not integrate across magnitude
or depth and so those values must be provided.
Args:
magnitude (float):
Earthquake magnitude.
hyp_depth (float):
Hypocentral depth (km).
ar (float):
Aspect ratio (L/W).
mechanism (Mechanism):
A ps2ff.constants Mechanism instance.
mag_scaling (MagScaling):
A ps2ff.constants MagScaling instance.
n_repi (int):
Number of log-spaced Repi points to compute conversions.
min_repi (float):
Minimum Repi to compute conversions (km).
max_repi (float):
Maximum Repi to compute conversions (km).
nxny (int):
Number of integration steps in the x/y direction for floating,
the rupture plane around the hypocenter.
n_theta (int):
Number of integration steps for theta. Default value of 9 is
an angular step size of 45 deg since it goes from 0 to 360 deg.
n_dip (int)
Number of integration steps for dip. Default value of 3 gives
a step size of 45 deg for the default range of 0 to 90 deg.
min_dip (float):
Minimum dip for integration (rad).
max_dip (float):
Maximum dip for integration (rad).
n_eps (int):
Number of integration steps for mag-area relationship.
trunc (float):
Truncation level for normal distribution of log area conditioned
on magnitude.
Retunrs:
tuple: All arrays of length n_repi:
- Repi (km).
- Average Rjb conditioned on Repi, M, and Zhyp.
- Average Rrup conditioned on Repi, M, and Zhyp.
- Variance of Rjb conditioned on Repi, M, and Zhyp.
- Variance of Rrup conditioned on Repi, M, and Zhyp.
"""
start_time = time.time()

# Check that mag_scaling and mechanism are the
# correct type
if not isinstance(mag_scaling, MagScaling):
raise TypeError(
'mag_scaling must bee a MagScaling instance')
if not isinstance(mechanism, Mechanism):
raise TypeError(
'mechanism must bee a Mechanism instance')

# Create conf dictionary.
conf = {}
conf['M'] = magnitude
conf['zhyp'] = hyp_depth
conf['AR'] = ar
conf['mech'] = mechanism
conf['rup_dim_model'] = mag_scaling

conf['nxny'] = nxny

conf['ndip'] = n_dip
conf['mindip'] = min_dip
conf['maxdip'] = max_dip

conf['neps'] = n_eps
conf['trunc'] = trunc

conf['min_seis_depth'] = 0.0
conf['max_seis_depth'] = 35.0 # not used?
conf['bytheta'] = False

repi_min = 0.1
repi_max = 1000
repi = np.logspace(
np.log10(repi_min), np.log10(repi_max), n_repi)

Rrup_var = np.zeros_like(repi)
Rrup_avg = np.zeros_like(repi)
Rjb_var = np.zeros_like(repi)
Rjb_avg = np.zeros_like(repi)

for i in range(n_repi):
Rrup_var[i], Rrup_avg[i], Rjb_var[i], Rjb_avg[i] = \
single_event_inner_loop(conf, repi[i], ntheta=n_theta)

print("Total execution time %f seconds." % (time.time() - start_time))

return repi, Rjb_avg, Rrup_avg, Rjb_var, Rrup_var
96 changes: 48 additions & 48 deletions tests/RjbRrup_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,18 +7,18 @@
def test_rjb_WC94():
cmd = "bin/run_ps2ff -w rjb tests/config/test_WC94.ini"
rc, so, se = get_command_output(cmd)
r1 = pd.DataFrame.from_csv(
"tests/data/test_Rjb_WC94_mechA_ar1p7_seis0_20_Ratios.csv",
header=6)
v1 = pd.DataFrame.from_csv(
"tests/data/test_Rjb_WC94_mechA_ar1p7_seis0_20_Var.csv",
header=6)
r2 = pd.DataFrame.from_csv(
"TestData/Rjb_WC94_mechA_ar1p7_seis0_20_Ratios.csv",
header=6)
v2 = pd.DataFrame.from_csv(
"TestData/Rjb_WC94_mechA_ar1p7_seis0_20_Var.csv",
header=6)
r1 = pd.read_csv(
"tests/data/test_Rjb_WC94_mechA_ar1p7_seis0_20_Ratios.csv",
header=6)
v1 = pd.read_csv(
"tests/data/test_Rjb_WC94_mechA_ar1p7_seis0_20_Var.csv",
header=6)
r2 = pd.read_csv(
"TestData/Rjb_WC94_mechA_ar1p7_seis0_20_Ratios.csv",
header=6)
v2 = pd.read_csv(
"TestData/Rjb_WC94_mechA_ar1p7_seis0_20_Var.csv",
header=6)

pd.util.testing.assert_frame_equal(r1, r2)
pd.util.testing.assert_frame_equal(v1, v2)
Expand All @@ -30,18 +30,18 @@ def test_rjb_WC94():
def test_rjb_S14():
cmd = "bin/run_ps2ff -w rjb tests/config/test_S14.ini"
rc, so, se = get_command_output(cmd)
r1 = pd.DataFrame.from_csv(
"tests/data/test_Rjb_S14_mechA_ar1p7_seis0_20_Ratios.csv",
header=6)
v1 = pd.DataFrame.from_csv(
"tests/data/test_Rjb_S14_mechA_ar1p7_seis0_20_Var.csv",
header=6)
r2 = pd.DataFrame.from_csv(
"TestData/Rjb_S14_mechA_ar1p7_seis0_20_Ratios.csv",
header=6)
v2 = pd.DataFrame.from_csv(
"TestData/Rjb_S14_mechA_ar1p7_seis0_20_Var.csv",
header=6)
r1 = pd.read_csv(
"tests/data/test_Rjb_S14_mechA_ar1p7_seis0_20_Ratios.csv",
header=6)
v1 = pd.read_csv(
"tests/data/test_Rjb_S14_mechA_ar1p7_seis0_20_Var.csv",
header=6)
r2 = pd.read_csv(
"TestData/Rjb_S14_mechA_ar1p7_seis0_20_Ratios.csv",
header=6)
v2 = pd.read_csv(
"TestData/Rjb_S14_mechA_ar1p7_seis0_20_Var.csv",
header=6)

pd.util.testing.assert_frame_equal(r1, r2)
pd.util.testing.assert_frame_equal(v1, v2)
Expand All @@ -53,18 +53,18 @@ def test_rjb_S14():
def test_rrup_S14():
cmd = "bin/run_ps2ff -w rrup tests/config/test_2_S14.ini"
rc, so, se = get_command_output(cmd)
r1 = pd.DataFrame.from_csv(
"tests/data/test_Rrup_S14_mechA_ar2p0_seis0_15_Ratios.csv",
header=6)
v1 = pd.DataFrame.from_csv(
"tests/data/test_Rrup_S14_mechA_ar2p0_seis0_15_Var.csv",
header=6)
r2 = pd.DataFrame.from_csv(
"TestData/Rrup_S14_mechA_ar2p0_seis0_15_Ratios.csv",
header=6)
v2 = pd.DataFrame.from_csv(
"TestData/Rrup_S14_mechA_ar2p0_seis0_15_Var.csv",
header=6)
r1 = pd.read_csv(
"tests/data/test_Rrup_S14_mechA_ar2p0_seis0_15_Ratios.csv",
header=6)
v1 = pd.read_csv(
"tests/data/test_Rrup_S14_mechA_ar2p0_seis0_15_Var.csv",
header=6)
r2 = pd.read_csv(
"TestData/Rrup_S14_mechA_ar2p0_seis0_15_Ratios.csv",
header=6)
v2 = pd.read_csv(
"TestData/Rrup_S14_mechA_ar2p0_seis0_15_Var.csv",
header=6)

pd.util.testing.assert_frame_equal(r1, r2)
pd.util.testing.assert_frame_equal(v1, v2)
Expand All @@ -76,18 +76,18 @@ def test_rrup_S14():
def test_rrup_WC94():
cmd = "bin/run_ps2ff -w rrup tests/config/test_2_WC94.ini"
rc, so, se = get_command_output(cmd)
r1 = pd.DataFrame.from_csv(
"tests/data/test_Rrup_WC94_mechA_ar2p0_seis0_15_Ratios.csv",
header=6)
v1 = pd.DataFrame.from_csv(
"tests/data/test_Rrup_WC94_mechA_ar2p0_seis0_15_Var.csv",
header=6)
r2 = pd.DataFrame.from_csv(
"TestData/Rrup_WC94_mechA_ar2p0_seis0_15_Ratios.csv",
header=6)
v2 = pd.DataFrame.from_csv(
"TestData/Rrup_WC94_mechA_ar2p0_seis0_15_Var.csv",
header=6)
r1 = pd.read_csv(
"tests/data/test_Rrup_WC94_mechA_ar2p0_seis0_15_Ratios.csv",
header=6)
v1 = pd.read_csv(
"tests/data/test_Rrup_WC94_mechA_ar2p0_seis0_15_Var.csv",
header=6)
r2 = pd.read_csv(
"TestData/Rrup_WC94_mechA_ar2p0_seis0_15_Ratios.csv",
header=6)
v2 = pd.read_csv(
"TestData/Rrup_WC94_mechA_ar2p0_seis0_15_Var.csv",
header=6)

pd.util.testing.assert_frame_equal(r1, r2)
pd.util.testing.assert_frame_equal(v1, v2)
Expand Down
Loading

0 comments on commit 19ebf28

Please sign in to comment.