Skip to content
GitLab
Menu
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Shakemap
shakyground2
Commits
90456796
Commit
90456796
authored
Mar 04, 2021
by
Marius Kriegerowski
Browse files
linting
parent
cb86d932
Pipeline
#20116
passed with stage
in 6 minutes and 8 seconds
Changes
11
Pipelines
2
Hide whitespace changes
Inline
Side-by-side
Makefile
View file @
90456796
SOURCES
=
shakyground2
*
.py
SOURCES
=
shakyground2
tests
*
.py
LENGTH
=
96
check
:
$(SOURCES)
...
...
shakyground2/earthquake.py
View file @
90456796
...
...
@@ -204,7 +204,7 @@ class Earthquake(object):
width
:
float
,
lsd
:
float
,
usd
:
float
=
0.0
,
hypo_loc
:
Tuple
[
float
,
float
]
=
(
0.5
,
0.5
)
hypo_loc
:
Tuple
[
float
,
float
]
=
(
0.5
,
0.5
)
,
):
"""
From a set of rupture properties returns a planar surface whose dimensions
...
...
@@ -266,12 +266,8 @@ class Earthquake(object):
updip_dir
=
(
dip
-
90.0
)
%
360
mid_left
=
centroid
.
point_at
(
left_length
,
0.0
,
(
strike
+
180.0
)
%
360.0
)
mid_right
=
centroid
.
point_at
(
right_length
,
0.0
,
strike
)
top_left
=
mid_left
.
point_at
(
updip_surface_length
,
-
updip_depth_change
,
updip_dir
)
top_right
=
mid_right
.
point_at
(
updip_surface_length
,
-
updip_depth_change
,
updip_dir
)
top_left
=
mid_left
.
point_at
(
updip_surface_length
,
-
updip_depth_change
,
updip_dir
)
top_right
=
mid_right
.
point_at
(
updip_surface_length
,
-
updip_depth_change
,
updip_dir
)
bottom_left
=
mid_left
.
point_at
(
downdip_surface_length
,
downdip_depth_change
,
downdip_dir
)
...
...
@@ -288,7 +284,7 @@ class Earthquake(object):
extended_error_message
=
[
"Rupture surface failed to build with the following properties:"
,
"Strike: %.2f, Dip: %.2f, Length: %.2f, Width: %.2f, Hypo. Pos: %s"
%
(
strike
,
dip
,
length
,
width
,
str
(
hypo_loc
))
%
(
strike
,
dip
,
length
,
width
,
str
(
hypo_loc
))
,
]
for
pnt
in
[
top_left
,
top_right
,
bottom_right
,
bottom_left
]:
extended_error_message
.
append
(
str
(
pnt
))
...
...
shakyground2/magnitude_scaling_relations.py
View file @
90456796
...
...
@@ -12,10 +12,17 @@ class BaseScalingRelation(metaclass=abc.ABCMeta):
Abstract base class representing an implementation of a magnitude scaling relation
"""
@
abc
.
abstractmethod
def
get_rupture_dimensions
(
self
,
magnitude
:
float
,
rake
:
float
=
0.0
,
dip
:
float
=
90.0
,
aspect
:
float
=
1.0
,
thickness
:
float
=
20.0
,
epsilon
:
float
=
0.0
)
->
Tuple
[
float
,
float
,
float
]:
def
get_rupture_dimensions
(
self
,
magnitude
:
float
,
rake
:
float
=
0.0
,
dip
:
float
=
90.0
,
aspect
:
float
=
1.0
,
thickness
:
float
=
20.0
,
epsilon
:
float
=
0.0
,
)
->
Tuple
[
float
,
float
,
float
]:
"""
Args:
magnitude:
...
...
@@ -41,9 +48,15 @@ class PEERScalingRelation(BaseScalingRelation):
updated by the synthetic rupture generator
"""
def
get_rupture_dimensions
(
self
,
magnitude
:
float
,
rake
:
float
=
0.0
,
dip
:
float
=
90.0
,
aspect
:
float
=
1.0
,
thickness
:
float
=
20.0
,
epsilon
:
float
=
0.0
)
->
Tuple
[
float
,
float
,
float
]:
def
get_rupture_dimensions
(
self
,
magnitude
:
float
,
rake
:
float
=
0.0
,
dip
:
float
=
90.0
,
aspect
:
float
=
1.0
,
thickness
:
float
=
20.0
,
epsilon
:
float
=
0.0
,
)
->
Tuple
[
float
,
float
,
float
]:
"""
Returns the area (km^2), length (km) and width (km) of the rupture using the PEER
Magnitude-Area scaling relation and constrained by the crustal thickness
...
...
@@ -80,6 +93,7 @@ class StrasserEtAl2010Interface(PEERScalingRelation):
Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude", Seismological
Research Letters, 81: 941 - 950, doi:10.1785/gssrl.81.6.941
"""
@
staticmethod
def
get_area
(
magnitude
:
float
)
->
float
:
"""
...
...
@@ -97,6 +111,7 @@ class StrasserEtAl2010Inslab(PEERScalingRelation):
Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude", Seismological
Research Letters, 81: 941 - 950, doi:10.1785/gssrl.81.6.941
"""
@
staticmethod
def
get_area
(
magnitude
:
float
)
->
float
:
"""
...
...
@@ -117,27 +132,39 @@ class Stafford2014(BaseScalingRelation):
# Model coefficients for the specific style of faulting
COEFFS
=
{
# Coefficients from Table 1
1
:
{
"U"
:
{
"a0"
:
-
27.4922
,
"a1"
:
4.6656
,
"a2"
:
-
0.2033
},
1
:
{
"U"
:
{
"a0"
:
-
27.4922
,
"a1"
:
4.6656
,
"a2"
:
-
0.2033
},
"SS"
:
{
"a0"
:
-
30.8395
,
"a1"
:
5.4184
,
"a2"
:
-
0.3044
},
"N"
:
{
"a0"
:
-
36.9770
,
"a1"
:
6.3070
,
"a2"
:
-
0.1696
},
"R"
:
{
"a0"
:
-
35.8239
,
"a1"
:
5.0680
,
"a2"
:
-
0.0457
}},
"R"
:
{
"a0"
:
-
35.8239
,
"a1"
:
5.0680
,
"a2"
:
-
0.0457
},
},
# Coefficients from Table 2
2
:
{
"U"
:
{
"b0"
:
-
2.300
,
"b1"
:
0.7167
,
"sigma"
:
0.2337
},
2
:
{
"U"
:
{
"b0"
:
-
2.300
,
"b1"
:
0.7167
,
"sigma"
:
0.2337
},
"SS"
:
{
"b0"
:
-
2.300
,
"b1"
:
0.7167
,
"sigma"
:
0.2337
},
"N"
:
{
"b0"
:
-
4.1055
,
"b1"
:
1.0370
,
"sigma"
:
0.2509
},
"R"
:
{
"b0"
:
-
3.8300
,
"b1"
:
0.9982
,
"sigma"
:
0.2285
}},
"R"
:
{
"b0"
:
-
3.8300
,
"b1"
:
0.9982
,
"sigma"
:
0.2285
},
},
# Coefficients from Table 3
3
:
{
"U"
:
{
"gamma0"
:
-
9.3137
,
"sigma"
:
0.3138
,
"rho"
:
0.3104
},
3
:
{
"U"
:
{
"gamma0"
:
-
9.3137
,
"sigma"
:
0.3138
,
"rho"
:
0.3104
},
"SS"
:
{
"gamma0"
:
-
9.3137
,
"sigma"
:
0.3138
,
"rho"
:
0.3104
},
"N"
:
{
"gamma0"
:
-
9.2483
,
"sigma"
:
0.3454
,
"rho"
:
0.4336
},
"R"
:
{
"gamma0"
:
-
9.2749
,
"sigma"
:
0.2534
,
"rho"
:
0.1376
}},
"R"
:
{
"gamma0"
:
-
9.2749
,
"sigma"
:
0.2534
,
"rho"
:
0.1376
},
},
# Coefficients from Table 4
4
:
{
"U"
:
0.7574
,
"SS"
:
0.7574
,
"N"
:
0.8490
,
"R"
:
0.7496
}
4
:
{
"U"
:
0.7574
,
"SS"
:
0.7574
,
"N"
:
0.8490
,
"R"
:
0.7496
}
,
}
def
get_rupture_dimensions
(
self
,
magnitude
:
float
,
rake
:
float
=
0.0
,
dip
:
float
=
90.0
,
aspect
:
float
=
1.0
,
thickness
:
float
=
20.0
,
epsilon
:
float
=
0.0
)
->
Tuple
[
float
,
float
,
float
]:
def
get_rupture_dimensions
(
self
,
magnitude
:
float
,
rake
:
float
=
0.0
,
dip
:
float
=
90.0
,
aspect
:
float
=
1.0
,
thickness
:
float
=
20.0
,
epsilon
:
float
=
0.0
,
)
->
Tuple
[
float
,
float
,
float
]:
"""
Gets the rupture dimensions from for the given magnitude subject to the physical
constaints.
...
...
@@ -193,8 +220,9 @@ class Stafford2014(BaseScalingRelation):
# normal
return
"N"
def
get_rupture_width
(
self
,
magnitude
:
float
,
dip
:
float
,
sof
:
str
,
thickness
:
float
)
->
\
Tuple
[
float
,
float
,
float
,
float
]:
def
get_rupture_width
(
self
,
magnitude
:
float
,
dip
:
float
,
sof
:
str
,
thickness
:
float
)
->
Tuple
[
float
,
float
,
float
,
float
]:
"""
Returns the parameters needed to define the censored rupture width
distribution defined in equations 8 to 14
...
...
@@ -206,44 +234,45 @@ class Stafford2014(BaseScalingRelation):
"""
# Gets the probability of a full width rupture
rw_max
=
thickness
/
sin
(
radians
(
dip
))
z_i
=
self
.
COEFFS
[
1
][
sof
][
"a0"
]
+
\
self
.
COEFFS
[
1
][
sof
][
"a1"
]
*
magnitude
+
\
self
.
COEFFS
[
1
][
sof
][
"a2"
]
*
rw_max
z_i
=
(
self
.
COEFFS
[
1
][
sof
][
"a0"
]
+
self
.
COEFFS
[
1
][
sof
][
"a1"
]
*
magnitude
+
self
.
COEFFS
[
1
][
sof
][
"a2"
]
*
rw_max
)
# Probability of rupturing full seismogenic thickness from logistic regression
p_i
=
1.0
/
(
1.0
+
exp
(
-
z_i
))
# Median width from equation 16
ln_rw
=
self
.
COEFFS
[
2
][
sof
][
"b0"
]
+
\
self
.
COEFFS
[
2
][
sof
][
"b1"
]
*
magnitude
ln_rw
=
self
.
COEFFS
[
2
][
sof
][
"b0"
]
+
self
.
COEFFS
[
2
][
sof
][
"b1"
]
*
magnitude
# Equations 18 - 20
phi_rw
=
(
log
(
rw_max
)
-
ln_rw
)
/
self
.
COEFFS
[
2
][
sof
][
"sigma"
]
phi_rw_ncdf
=
norm
.
cdf
(
phi_rw
)
ln_rw_trunc
=
ln_rw
-
self
.
COEFFS
[
2
][
sof
][
"sigma"
]
*
\
(
norm
.
pdf
(
phi_rw
)
/
phi_rw_ncdf
)
ln_rw_trunc
=
ln_rw
-
self
.
COEFFS
[
2
][
sof
][
"sigma"
]
*
(
norm
.
pdf
(
phi_rw
)
/
phi_rw_ncdf
)
mean_rw
=
p_i
*
log
(
rw_max
)
+
(
1.0
-
p_i
)
*
ln_rw_trunc
# Equations 21 - 22
stddev_rw
=
self
.
_get_rupture_width_sigma
(
self
.
COEFFS
[
2
][
sof
][
"sigma"
],
phi_rw
,
phi_rw_ncdf
,
p_i
)
stddev_rw
=
self
.
_get_rupture_width_sigma
(
self
.
COEFFS
[
2
][
sof
][
"sigma"
],
phi_rw
,
phi_rw_ncdf
,
p_i
)
return
mean_rw
,
stddev_rw
,
rw_max
,
p_i
@
staticmethod
def
_get_rupture_width_sigma
(
sigma
:
float
,
phi_rw
:
float
,
phi_rw_ncdf
:
float
,
p_i
:
float
)
->
float
:
def
_get_rupture_width_sigma
(
sigma
:
float
,
phi_rw
:
float
,
phi_rw_ncdf
:
float
,
p_i
:
float
)
->
float
:
"""
Returns the variabiliy in the rupture width described by equations 21 and 22
"""
denom
=
sqrt
(
2.0
*
pi
)
*
phi_rw_ncdf
if
phi_rw_ncdf
>=
0.0
:
elem1
=
sqrt
(
pi
/
2.
)
*
(
1.0
+
chi2
.
cdf
(
phi_rw
,
3
))
elem1
=
sqrt
(
pi
/
2.
0
)
*
(
1.0
+
chi2
.
cdf
(
phi_rw
,
3
))
else
:
elem1
=
sqrt
(
pi
/
2.
)
*
(
1.0
-
chi2
.
cdf
(
phi_rw
,
3
))
elem1
=
sqrt
(
pi
/
2.
0
)
*
(
1.0
-
chi2
.
cdf
(
phi_rw
,
3
))
elem2
=
exp
(
-
(
phi_rw
**
2
))
/
denom
sigma_truncated
=
sqrt
(((
sigma
**
2.
)
/
denom
)
*
(
elem1
-
elem2
))
sigma_truncated
=
sqrt
(((
sigma
**
2.
0
)
/
denom
)
*
(
elem1
-
elem2
))
return
(
1.0
-
p_i
)
*
sigma_truncated
def
get_rupture_area
(
self
,
magnitude
:
float
,
sof
:
str
,
rw_max
:
float
,
sigma_lnrw
:
float
)
->
Tuple
[
float
,
float
]:
def
get_rupture_area
(
self
,
magnitude
:
float
,
sof
:
str
,
rw_max
:
float
,
sigma_lnrw
:
float
)
->
Tuple
[
float
,
float
]:
"""
Returns the rupture area conditioned upon the maximum rupture width and style of
faulting
...
...
@@ -258,15 +287,16 @@ class Stafford2014(BaseScalingRelation):
ln_ra: Natural logarithm of rupture area
sigma_ra: Standard deviation of the natural logarithm of rupture area
"""
mw_crit
=
(
log
(
rw_max
)
-
self
.
COEFFS
[
2
][
sof
][
"b0"
])
/
\
self
.
COEFFS
[
2
][
sof
][
"b1"
]
ln_ra
=
self
.
COEFFS
[
3
][
sof
][
"gamma0"
]
+
log
(
10.
)
*
magnitude
mw_crit
=
(
log
(
rw_max
)
-
self
.
COEFFS
[
2
][
sof
][
"b0"
])
/
self
.
COEFFS
[
2
][
sof
][
"b1"
]
ln_ra
=
self
.
COEFFS
[
3
][
sof
][
"gamma0"
]
+
log
(
10.0
)
*
magnitude
if
magnitude
>
mw_crit
:
# Equation 23
ln_ra
-=
(
(
log
(
10.
)
/
4.
)
*
(
magnitude
-
mw_crit
)
)
ln_ra
-=
(
log
(
10.
0
)
/
4.
0
)
*
(
magnitude
-
mw_crit
)
# Get the sigma log rupture area (equation 28)
sigma
=
self
.
COEFFS
[
3
][
sof
][
"sigma"
]
sigma_ra
=
sqrt
(
(
sigma
**
2.
)
+
(
sigma_lnrw
**
2.
)
+
(
2.0
*
self
.
COEFFS
[
3
][
sof
][
"rho"
]
*
sigma_lnrw
*
sigma
))
(
sigma
**
2.0
)
+
(
sigma_lnrw
**
2.0
)
+
(
2.0
*
self
.
COEFFS
[
3
][
sof
][
"rho"
]
*
sigma_lnrw
*
sigma
)
)
return
ln_ra
,
sigma_ra
shakyground2/rupture_mechanism.py
View file @
90456796
...
...
@@ -17,6 +17,7 @@ class DIP_RANGES(Enum):
for each style-of-faulting is provided and they are assumed to be uniformly
distributed.
"""
R
=
(
20.0
,
45.0
)
SS
=
(
75.0
,
91.0
)
N
=
(
45.0
,
75.0
)
...
...
@@ -40,6 +41,7 @@ class RuptureMechanism(object):
of :class:`openquake.hazardlib.pmf.PMF`
"""
def
__init__
(
self
,
mechanism
:
Optional
[
PMF
]
=
None
):
if
mechanism
:
assert
isinstance
(
mechanism
,
PMF
)
...
...
@@ -76,8 +78,12 @@ class RuptureMechanism(object):
return
strikes
,
dips
,
rakes
@
classmethod
def
from_strike_dip_rake
(
cls
,
strike
:
Optional
[
float
]
=
None
,
dip
:
Optional
[
float
]
=
None
,
rake
:
Optional
[
float
]
=
None
)
->
RuptureMechanism
:
def
from_strike_dip_rake
(
cls
,
strike
:
Optional
[
float
]
=
None
,
dip
:
Optional
[
float
]
=
None
,
rake
:
Optional
[
float
]
=
None
,
)
->
RuptureMechanism
:
"""
Constructs the rupture mechanism distribution from a simple strike, dip
and rake combination, permitting None values if undefined
...
...
@@ -87,13 +93,14 @@ class RuptureMechanism(object):
dip: Dip of fault (in decimal degrees between 0 and 90)
rake: Rake of fault (in decimal degrees between -180 and 180)
"""
return
cls
(
cls
.
build_mechanism_distribution
(
valid
.
strike
(
strike
),
valid
.
dip
(
dip
),
valid
.
rake
(
rake
)))
return
cls
(
cls
.
build_mechanism_distribution
(
valid
.
strike
(
strike
),
valid
.
dip
(
dip
),
valid
.
rake
(
rake
)
)
)
@
classmethod
def
from_nodal_planes
(
cls
,
nodal_plane_1
:
Dict
,
nodal_plane_2
:
Dict
)
\
->
RuptureMechanism
:
def
from_nodal_planes
(
cls
,
nodal_plane_1
:
Dict
,
nodal_plane_2
:
Dict
)
->
RuptureMechanism
:
"""
Constructs the rupture mechanism distribution from either one single
valid rupture plane or two valud nodal planes
...
...
@@ -108,13 +115,21 @@ class RuptureMechanism(object):
strike_2
=
valid
.
strike
(
nodal_plane_2
[
"strike"
])
dip_2
=
valid
.
dip
(
nodal_plane_2
[
"dip"
])
rake_2
=
valid
.
rake
(
nodal_plane_2
[
"rake"
])
return
cls
(
PMF
([(
0.5
,
NodalPlane
(
strike_1
,
dip_1
,
rake_1
)),
(
0.5
,
NodalPlane
(
strike_2
,
dip_2
,
rake_2
))]))
return
cls
(
PMF
(
[
(
0.5
,
NodalPlane
(
strike_1
,
dip_1
,
rake_1
)),
(
0.5
,
NodalPlane
(
strike_2
,
dip_2
,
rake_2
)),
]
)
)
@
staticmethod
def
build_mechanism_distribution
(
strike
:
Optional
[
float
]
=
None
,
dip
:
Optional
[
float
]
=
None
,
rake
:
Optional
[
float
]
=
None
)
->
PMF
:
def
build_mechanism_distribution
(
strike
:
Optional
[
float
]
=
None
,
dip
:
Optional
[
float
]
=
None
,
rake
:
Optional
[
float
]
=
None
,
)
->
PMF
:
"""
Builds a mechanism distribution from a partial (or complete) characterisation of the
rupture mechanism
...
...
shakyground2/site_model.py
View file @
90456796
...
...
@@ -23,7 +23,7 @@ SITE_PROPERTIES = {
"xvf"
:
np
.
float64
,
# Distance (km) to the volcanic front
"backarc"
:
np
.
bool
,
# Site is in the subduction backarc (True) or else
"region"
:
np
.
int32
,
# Region to which the site belongs
"geology"
:
(
np
.
string_
,
20
)
# Geological classification for the site
"geology"
:
(
np
.
string_
,
20
)
,
# Geological classification for the site
}
...
...
@@ -218,15 +218,11 @@ class SiteModel(object):
if
key
not
in
dframe
.
columns
:
# Use the default
if
SITE_PROPERTIES
[
key
]
in
((
np
.
bytes_
,
20
),):
site_array
[
key
]
=
cls
.
_get_string_array
(
nsites
,
value
,
SITE_PROPERTIES
[
key
]
)
site_array
[
key
]
=
cls
.
_get_string_array
(
nsites
,
value
,
SITE_PROPERTIES
[
key
])
elif
SITE_PROPERTIES
[
key
]
==
np
.
bool
:
site_array
[
key
]
=
cls
.
_get_boolean_array
(
nsites
,
value
)
else
:
site_array
[
key
]
=
value
*
np
.
ones
(
nsites
,
dtype
=
SITE_PROPERTIES
[
key
]
)
site_array
[
key
]
=
value
*
np
.
ones
(
nsites
,
dtype
=
SITE_PROPERTIES
[
key
])
return
cls
(
site_array
)
@
staticmethod
...
...
@@ -277,9 +273,7 @@ class SiteModel(object):
if
japan
:
c1
=
412.0
**
2.0
c2
=
1360.0
**
2.0
return
np
.
exp
(
(
-
5.23
/
2.0
)
*
np
.
log
((
np
.
power
(
vs30
,
2.0
)
+
c1
)
/
(
c2
+
c1
))
)
return
np
.
exp
((
-
5.23
/
2.0
)
*
np
.
log
((
np
.
power
(
vs30
,
2.0
)
+
c1
)
/
(
c2
+
c1
)))
else
:
c1
=
571.0
**
4.0
c2
=
1360.0
**
4.0
...
...
shakyground2/valid.py
View file @
90456796
...
...
@@ -4,8 +4,7 @@ consistent quantities
"""
from
typing
import
Tuple
,
Optional
,
Union
from
openquake.hazardlib.geo.nodalplane
import
NodalPlane
from
shakyground2.magnitude_scaling_relations
import
(
BaseScalingRelation
,
PEERScalingRelation
)
from
shakyground2.magnitude_scaling_relations
import
BaseScalingRelation
,
PEERScalingRelation
def
longitude
(
lon
:
float
)
->
float
:
...
...
@@ -93,8 +92,10 @@ def seismogenic_thickness(
usd
=
positive_float
(
upper_seismo_depth
,
"upper seismogenic depth"
)
lsd
=
positive_float
(
lower_seismo_depth
,
"lower seismogenic depth"
)
if
lsd
<
usd
:
raise
ValueError
(
"Lower seismogenic depth %.2f km shallower than upper seismogenic "
"depth %.2f km"
%
(
lsd
,
usd
))
raise
ValueError
(
"Lower seismogenic depth %.2f km shallower than upper seismogenic "
"depth %.2f km"
%
(
lsd
,
usd
)
)
return
usd
,
lsd
...
...
@@ -105,8 +106,9 @@ def hypocenter_position(hypo_pos: Tuple[float, float]) -> Tuple[float, float]:
"""
along_strike
,
down_dip
=
hypo_pos
if
along_strike
<
0.0
or
along_strike
>
1.0
:
raise
ValueError
(
"Along strike position %.3f should be in the range 0 to 1"
%
along_strike
)
raise
ValueError
(
"Along strike position %.3f should be in the range 0 to 1"
%
along_strike
)
if
down_dip
<
0.0
or
down_dip
>
1.0
:
raise
ValueError
(
"Down dip position %.3f should be in the range 0 to 1"
%
down_dip
)
return
along_strike
,
down_dip
...
...
@@ -121,6 +123,7 @@ def scaling_relation(msr: Optional[BaseScalingRelation]):
# If no relation is defined then use the default
return
PEERScalingRelation
()
if
not
isinstance
(
msr
,
BaseScalingRelation
):
raise
TypeError
(
"Magnitude Scaling Relation %s not instance of BaseScalingRelation"
%
str
(
msr
))
raise
TypeError
(
"Magnitude Scaling Relation %s not instance of BaseScalingRelation"
%
str
(
msr
)
)
return
msr
tests/earthquake_test.py
View file @
90456796
...
...
@@ -11,7 +11,6 @@ from shakyground2.magnitude_scaling_relations import PEERScalingRelation
class
EarthquakeTestCase
(
unittest
.
TestCase
):
def
setUp
(
self
):
self
.
lon
=
35.0
self
.
lat
=
40.0
...
...
@@ -27,10 +26,18 @@ class EarthquakeTestCase(unittest.TestCase):
self
.
assertEqual
(
str
(
eq0
),
"XYZ (35.00000E, 40.00000N, 10.00 km) M 6.00"
)
# Minimal instantiation - no datetime
eq0
=
Earthquake
(
"XYZ"
,
self
.
lon
,
self
.
lat
,
self
.
depth
,
self
.
mag
,
date
=
date
(
2010
,
10
,
6
),
time
=
time
(
11
,
35
,
45
))
self
.
assertEqual
(
str
(
eq0
),
"XYZ 2010-10-06 11:35:45 (35.00000E, 40.00000N, 10.00 km) M 6.00"
)
eq0
=
Earthquake
(
"XYZ"
,
self
.
lon
,
self
.
lat
,
self
.
depth
,
self
.
mag
,
date
=
date
(
2010
,
10
,
6
),
time
=
time
(
11
,
35
,
45
),
)
self
.
assertEqual
(
str
(
eq0
),
"XYZ 2010-10-06 11:35:45 (35.00000E, 40.00000N, 10.00 km) M 6.00"
)
def
test_create_planar_surface_simple
(
self
):
# Test the creation of a simple planar surface that requires no
...
...
@@ -41,10 +48,10 @@ class EarthquakeTestCase(unittest.TestCase):
dip
=
90.0
usd
=
0.0
lsd
=
20.0
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
# Length should be approximately 10 km, but allow a slight precision difference
self
.
assertAlmostEqual
(
plane0
.
length
,
10.0
,
1
)
self
.
assertAlmostEqual
(
plane0
.
width
,
10.0
,
5
)
...
...
@@ -66,10 +73,10 @@ class EarthquakeTestCase(unittest.TestCase):
dip
=
90.0
usd
=
0.0
lsd
=
20.0
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
# Verify correct plane from top left and bottom right points
self
.
assertAlmostEqual
(
plane0
.
top_left
.
longitude
,
-
self
.
target_width
/
2.0
,
7
)
...
...
@@ -88,10 +95,10 @@ class EarthquakeTestCase(unittest.TestCase):
dip
=
90.0
usd
=
0.0
lsd
=
10.0
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
# Verify correct plane from top left and bottom right points
self
.
assertAlmostEqual
(
plane0
.
top_left
.
longitude
,
-
self
.
target_width
/
2.0
,
7
)
...
...
@@ -110,17 +117,16 @@ class EarthquakeTestCase(unittest.TestCase):
dip
=
90.0
usd
=
0.0
lsd
=
10.0
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
# Verify correct plane from top left and bottom right points
self
.
assertAlmostEqual
(
plane0
.
top_left
.
longitude
,
-
10.0
*
self
.
target_width
/
2.0
,
7
)
self
.
assertAlmostEqual
(
plane0
.
top_left
.
latitude
,
0.0
,
7
)
self
.
assertAlmostEqual
(
plane0
.
top_left
.
depth
,
0.0
,
7
)
self
.
assertAlmostEqual
(
plane0
.
bottom_right
.
longitude
,
10.0
*
self
.
target_width
/
2.0
,
7
)
self
.
assertAlmostEqual
(
plane0
.
bottom_right
.
longitude
,
10.0
*
self
.
target_width
/
2.0
,
7
)
self
.
assertAlmostEqual
(
plane0
.
bottom_right
.
latitude
,
0.0
,
7
)
self
.
assertAlmostEqual
(
plane0
.
bottom_right
.
depth
,
10.0
,
7
)
...
...
@@ -133,10 +139,10 @@ class EarthquakeTestCase(unittest.TestCase):
dip
=
60.0
usd
=
0.0
lsd
=
10.0
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
self
.
assertAlmostEqual
(
width
,
10.0
/
sin
(
radians
(
dip
)))
# Verify correct plane from top left and bottom right points - calculated by hand
...
...
@@ -154,34 +160,58 @@ class EarthquakeTestCase(unittest.TestCase):
strike
=
0.0
dip
=
90.0
# Create a simple surface
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
self
.
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
area
,
length
,
width
=
PEERScalingRelation
().
get_rupture_dimensions
(
self
.
mag
,
dip
=
dip
,
thickness
=
lsd
-
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
plane0
=
Earthquake
.
build_planar_surface
(
centroid
,
strike
,
dip
,
length
,
width
,
lsd
,
usd
)
# Instantiate earthquake object with known surface
eq0
=
Earthquake
(
"XYZ"
,
self
.
lon
,
self
.
lat
,
self
.
depth
,
self
.
mag
,
strike
=
strike
,
dip
=
dip
,
rake
=
0.0
,
upper_seismogenic_depth
=
usd
,
lower_seismogenic_depth
=
lsd
,
surface
=
plane0
)
eq0
=
Earthquake
(
"XYZ"
,
self
.
lon
,
self
.
lat
,
self
.
depth
,
self
.
mag
,
strike
=
strike
,
dip
=
dip
,
rake
=
0.0
,
upper_seismogenic_depth
=
usd
,
lower_seismogenic_depth
=
lsd
,
surface
=
plane0
,
)
# Generate rupture
self
.
assertAlmostEqual
(
eq0
.
rupture
.
mag
,
self
.
mag
)
self
.
assertAlmostEqual
(
eq0
.
rupture
.
rake
,
0.0
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_lons
,
eq0
.
rupture
.
surface
.
corner_lons
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_lats
,
eq0
.
rupture
.
surface
.
corner_lats
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_depths
,
eq0
.
rupture
.
surface
.
corner_depths
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_lons
,
eq0
.
rupture
.
surface
.
corner_lons
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_lats
,
eq0
.
rupture
.
surface
.
corner_lats
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_depths
,
eq0
.
rupture
.
surface
.
corner_depths
)
# Instantiate earthquake object without surface
eq1
=
Earthquake
(
"XYZ"
,
self
.
lon
,
self
.
lat
,
self
.
depth
,
self
.
mag
,
strike
=
strike
,
dip
=
dip
,
rake
=
0.0
,
upper_seismogenic_depth
=
usd
,
lower_seismogenic_depth
=
lsd
)
eq1
=
Earthquake
(
"XYZ"
,
self
.
lon
,
self
.
lat
,
self
.
depth
,
self
.
mag
,
strike
=
strike
,
dip
=
dip
,
rake
=
0.0
,
upper_seismogenic_depth
=
usd
,
lower_seismogenic_depth
=
lsd
,
)
self
.
assertAlmostEqual
(
eq1
.
rupture
.
mag
,
self
.
mag
)
self
.
assertAlmostEqual
(
eq1
.
rupture
.
rake
,
0.0
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_lons
,
eq1
.
rupture
.
surface
.
corner_lons
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_lats
,
eq1
.
rupture
.
surface
.
corner_lats
)
np
.
testing
.
assert_array_almost_equal
(
plane0
.
corner_depths
,