Skip to content
GitLab
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
cd06865f
Commit
cd06865f
authored
Mar 12, 2021
by
g-weatherill
Browse files
Merge remote-tracking branch 'upstream/master'
parents
dce96698
e61ff9d3
Changes
11
Hide whitespace changes
Inline
Side-by-side
shakyground2/earthquake.py
View file @
cd06865f
...
...
@@ -2,12 +2,13 @@
Contains the core classes to represent an earthquake
"""
import
datetime
from
typing
import
Tuple
from
typing
import
Tuple
,
Optional
,
List
from
math
import
radians
,
sin
,
tan
from
openquake.hazardlib.geo
import
Point
,
PlanarSurface
from
openquake.hazardlib.geo.surface.base
import
BaseSurface
from
openquake.hazardlib.source.rupture
import
ParametricProbabilisticRupture
from
shakyground2
import
valid
from
shakyground2.rupture_mechanism
import
RuptureMechanism
class
Earthquake
(
object
):
...
...
@@ -28,6 +29,9 @@ class Earthquake(object):
Earthquake latitude in decimal degrees (-90, to 90)
depth:
Hypocentral depth (in km)
hypocenter:
Representation of the hypocenter as an instance of
:class:`openquake.hazardlib.geo.Point`
mag:
Earthquake magnitude
strike:
...
...
@@ -63,14 +67,15 @@ class Earthquake(object):
lat
,
hypo_depth
,
magnitude
,
strike
=
0.0
,
dip
=
90.0
,
rake
=
0.0
,
strike
=
None
,
dip
=
None
,
rake
=
None
,
aspect
=
1.0
,
mag_scaling_relation
=
None
,
upper_seismogenic_depth
=
0.0
,
lower_seismogenic_depth
=
1000.0
,
surface
=
None
,
focal_mechanism
=
None
,
hypocenter_position
=
(
0.5
,
0.5
),
date
=
None
,
time
=
None
,
...
...
@@ -121,6 +126,7 @@ class Earthquake(object):
self
.
lon
=
valid
.
longitude
(
lon
)
self
.
lat
=
valid
.
latitude
(
lat
)
self
.
depth
=
valid
.
positive_float
(
hypo_depth
,
"hypocentre depth"
)
self
.
hypocenter
=
Point
(
self
.
lon
,
self
.
lat
,
self
.
depth
)
self
.
mag
=
magnitude
self
.
strike
=
valid
.
strike
(
strike
)
self
.
dip
=
valid
.
dip
(
dip
)
...
...
@@ -138,8 +144,41 @@ class Earthquake(object):
self
.
time
=
time
assert
isinstance
(
surface
,
BaseSurface
)
or
surface
is
None
self
.
surface
=
surface
if
self
.
surface
:
# Can calculate rupture dimensions from the surface
self
.
area
=
self
.
surface
.
get_area
()
self
.
width
=
self
.
surface
.
get_width
()
self
.
length
=
self
.
area
/
self
.
width
else
:
# Default rupture dimensions to none to none
self
.
area
=
self
.
width
=
self
.
length
=
None
self
.
hypo_pos
=
valid
.
hypocenter_position
(
hypocenter_position
)
self
.
_rupture
=
None
# Get a valid focal mechanism with two nodal planes
self
.
focal_mechanism
=
valid
.
focal_mechanism
(
focal_mechanism
)
self
.
mechanism
=
self
.
_get_mechanism
()
self
.
tectonic_region
=
None
def
_get_mechanism
(
self
):
"""
Defines the focal mechanism according to three different cases:
1. A unique mechanism is defined explicitly from the strike, dip and rake
2. A pair of equiprobable mechanisms is defined from the focal mechanism
3. No mechanism is defined, in which case a global distribution is assumed
Returns:
Mechanism distribution as an instance of the
:class:`shakyground2.rupture_mechanism.RuptureMechanism`
"""
if
(
self
.
strike
is
not
None
)
and
(
self
.
dip
is
not
None
)
and
(
self
.
rake
is
not
None
):
# Fixed and fully defined rupture mechanism
return
RuptureMechanism
.
from_strike_dip_rake
(
self
.
strike
,
self
.
dip
,
self
.
rake
)
elif
self
.
focal_mechanism
:
# Rupture mechanism defines from nodal planes
return
RuptureMechanism
.
from_focal_mechanism
(
self
.
focal_mechanism
)
else
:
# Global distribution
return
RuptureMechanism
()
def
__repr__
(
self
):
# Returns a summary string of the event, with datetime if specified
...
...
@@ -156,6 +195,34 @@ class Earthquake(object):
self
.
id
,
self
.
lon
,
self
.
lat
,
self
.
depth
,
self
.
mag
)
def
get_maximum_distance_bbox
(
self
,
max_distance
:
Optional
[
float
]
=
None
)
->
List
:
"""
Defines a bounding box around the event up to a maximum distance.
Args:
max_distance: Maximum horizontal and vertical distance from the epicentre for the
bounding box to be defined. If None then a default bounding box
size is determined that is based on magnitude between Mw 3.0 (100 km)
and Mw 8.0 (1000 km)
Returns:
Bounding box as a list of [llon, llat, ulon, ulat]
"""
if
not
max_distance
:
# Scale from 100 km (for Mw 3 or less) to 1000 km for Mw >= 8.0
if
self
.
mag
<=
3.0
:
max_distance
=
100.0
elif
self
.
mag
>=
8.0
:
max_distance
=
1000.0
else
:
# Interpolate
max_distance
=
100.0
+
(
self
.
mag
-
3.0
)
*
(
900.0
/
5.0
)
# Define the bounding box from the specified maximum distance
north
=
self
.
hypocenter
.
point_at
(
max_distance
,
0.0
,
0.0
)
south
=
self
.
hypocenter
.
point_at
(
max_distance
,
0.0
,
180.0
)
east
=
self
.
hypocenter
.
point_at
(
max_distance
,
0.0
,
90.0
)
west
=
self
.
hypocenter
.
point_at
(
max_distance
,
0.0
,
270.0
)
return
[
west
.
longitude
,
south
.
latitude
,
east
.
longitude
,
north
.
latitude
]
@
property
def
rupture
(
self
):
"""
...
...
@@ -172,27 +239,6 @@ class Earthquake(object):
self
.
mag
,
self
.
rake
,
None
,
centroid
,
self
.
surface
,
1.0
,
None
)
return
self
.
_rupture
# Rupture is not defined, so needs to be constructed. Use simple case for now
# TODO To be replaced by the synthetic rupture generator operation
thickness
=
self
.
lsd
-
self
.
usd
# Get rupture dimensions
area
,
length
,
width
=
self
.
mag_scale_rel
.
get_rupture_dimensions
(
self
.
mag
,
aspect
=
self
.
aspect
,
dip
=
self
.
dip
,
thickness
=
thickness
)
# Build the rupture surface
surface
=
self
.
build_planar_surface
(
centroid
,
self
.
strike
,
self
.
dip
,
length
,
width
,
self
.
lsd
,
self
.
usd
,
self
.
hypo_pos
,
)
self
.
_rupture
=
ParametricProbabilisticRupture
(
self
.
mag
,
self
.
rake
,
None
,
centroid
,
surface
,
1.0
,
None
)
return
self
.
_rupture
@
staticmethod
...
...
shakyground2/rupture_mechanism.py
View file @
cd06865f
...
...
@@ -4,7 +4,7 @@ Class to manage the distribution of rupture mechanisms
from
__future__
import
annotations
import
numpy
as
np
from
enum
import
Enum
from
typing
import
Optional
,
Dict
,
Tuple
from
typing
import
Optional
,
Dict
,
Tuple
,
List
from
openquake.hazardlib.geo
import
NodalPlane
from
openquake.hazardlib.pmf
import
PMF
from
shakyground2
import
valid
...
...
@@ -56,6 +56,10 @@ class RuptureMechanism(object):
for
prob
,
nodal_plane
in
self
.
mechanism
.
data
:
yield
prob
,
nodal_plane
def
__len__
(
self
):
# Lenth corresponds to the number of mechanisms in the distribution
return
len
(
self
.
mechanism
.
data
)
def
sample
(
self
,
nsamples
:
int
)
->
Tuple
[
np
.
ndarray
,
np
.
ndarray
,
np
.
ndarray
]:
"""
Samples the focal mechanism distribution to return "nsamples" strike,
...
...
@@ -100,30 +104,48 @@ class RuptureMechanism(object):
)
@
classmethod
def
from_
nodal_planes
(
cls
,
nodal_plane_1
:
Dict
,
nodal_plane_2
:
Dict
)
->
RuptureMechanism
:
def
from_
focal_mechanism
(
cls
,
focal_mechanism
:
Dict
)
->
RuptureMechanism
:
"""
Constructs the
rupture
mechanism
distribution from either o
ne s
ingle
valid rupture plane or two valud nodal planes
Constructs the
focal
mechanism
from an evenly pair of nodal pla
ne
s
s
uch as that
of a focal mechanism
Args:
nodal_plane_1: Dictionary containing strike, dip and rake of the first nodel plane
nodal_plane_1: Dictionary containing strike, dip and rake of the second nodel plane
focal_mechanism: Dictionary containing two :class:
`openquake.hazardlib.geo.NodalPlane` objects, each labeled as "nodal_plane_1"
and "nodal_plane_2"
"""
strike_1
=
valid
.
strike
(
nodal_plane_1
[
"strike"
])
dip_1
=
valid
.
dip
(
nodal_plane_1
[
"dip"
])
rake_1
=
valid
.
rake
(
nodal_plane_1
[
"rake"
])
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
,
rak
e_1
)
),
(
0.5
,
NodalPlane
(
strike_2
,
dip_2
,
rak
e_2
)
),
(
0.5
,
focal_mechanism
[
"nodal_plan
e_1
"
]
),
(
0.5
,
focal_mechanism
[
"nodal_plan
e_2
"
]
),
]
)
)
@
classmethod
def
from_nodal_planes
(
cls
,
nodal_planes
:
List
,
probabilities
:
List
)
->
RuptureMechanism
:
"""
Constructs the rupture mechanism distribution from a list of nodal planes and their
associated probabilities
Args:
nodal_planes: Set of nodal planes as a list of dictionaries, eac containing strike,
dip and rake
probabilities: List of probabilities of the nodal planes (must sum to 1)
"""
assert
len
(
nodal_planes
)
==
len
(
probabilities
),
"Number of nodal planes not equal to number of probabilities"
assert
np
.
isclose
(
sum
(
probabilities
),
1.0
),
"Probabilities do not sum to 1.0 (sum = %.6f)"
%
sum
(
probabilities
)
mechanism_distribution
=
[]
for
prob
,
npl
in
zip
(
probabilities
,
nodal_planes
):
mechanism
=
valid
.
mechanism
(
npl
[
"strike"
],
npl
[
"dip"
],
npl
[
"rake"
])
mechanism_distribution
.
append
((
prob
,
mechanism
))
return
cls
(
PMF
(
mechanism_distribution
))
@
staticmethod
def
build_mechanism_distribution
(
strike
:
Optional
[
float
]
=
None
,
...
...
shakyground2/shakemap.py
0 → 100644
View file @
cd06865f
"""
Implements the core shakemap class
"""
import
h5py
import
numpy
as
np
import
warnings
from
typing
import
Dict
,
Optional
,
Tuple
,
List
from
openquake.hazardlib
import
const
,
imt
from
openquake.hazardlib.contexts
import
ContextMaker
from
shakyground2
import
valid
from
shakyground2.earthquake
import
Earthquake
from
shakyground2.site_model
import
SiteModel
from
shakyground2.synthetic_rupture_generator
import
FiniteRuptureSampler
class
Shakemap
(
object
):
"""
The core class for calculating shakemaps. The general workflow contained within
takes the input earthquake and site model in order to build the source, path and site
context objects (i.e. containers that hold the data needed for the ground motion models).
The contexts objects are built upon instantiation of the shakemap class to avoid
recalculation if further shakemaps are needed.
The Shakemap class also caches the information into an hdf5 binary object, if desired by
the user. This means that the shakemap information for an earthquake ccould be retrieved at
a later time without the need for recalculation.
"""
def
__init__
(
self
,
earthquake
:
Earthquake
,
site_model
:
SiteModel
,
ground_motion_model
:
Dict
,
tectonic_region
:
str
,
cache_file
:
Optional
[
str
]
=
None
,
num_rupture_samples
:
int
=
100
,
rdim
:
float
=
0.0
,
synth_dist_weights
:
List
=
[
0.25
,
0.25
,
0.25
,
0.25
],
synthetic_rupture_max_site_distance
:
float
=
200.0
,
synthetic_rupture_site_spacing
:
float
=
0.05
,
):
"""
Args:
earthquake:
Earthquake as an instance of the :class:shakyground2.earthquake.Earthquake
site_model:
Target sites used for calculation of the ground motion values, as instances of
:class:shakyground2.site_model.SiteModel
ground_motion_model:
Set of ground motion models and their respective weights as a dictionary
{"GMMs": [...], "weights": [...]}
tectonic_region:
Tectonic region to which the earthquake is assigned
cache_file:
Path to the hdf5 file for caching the results. It does not need to exist
but if the hdf5 file exists and contains an instance of the same event
(according to the earthquake.id) then the user will be warned here and the
existing shakemaps for that earthquake ID subsequently overwitten if the
user persists to implement the shakemap.
num_rupture_samples:
In the likely case that no rupture surface is available for the earthquake,
one will be generated by the :class:`shakyground2.synthetic_rupture_generator.
FiniteRuptureSampler`. This controls the number of samples to be used
rdim:
Exponent of the distance-dependent site weighing adopted by the :class:
`shakyground2.synthetic_rupture_generator.FiniteRuptureSampler`
synth_dist_weights: Can adjust the weight assigned to the different distance
metrics used by the FiniteRuptureSampler
synthetic_rupture max_site_distance: In the case that the target sites for the
FiniteRuptureSampler need to be built from the
bounding box, this defines the maximum
distances of the target sites to define the
bounding box
synthetic_rupture_site_spacing: In the case that the target sites for the
FiniteRuptureSampler need to be built from the
bounding box, this defines the site spacing of the
sites
"""
self
.
earthquake
=
earthquake
self
.
site_model
=
site_model
self
.
tectonic_region
=
tectonic_region
self
.
ground_motion_model
=
ground_motion_model
self
.
num_rupture_samples
=
num_rupture_samples
self
.
rdim
=
rdim
self
.
synthetic_rupture_distance_weights
=
synth_dist_weights
if
cache_file
:
self
.
caching
=
True
self
.
cache_file
=
cache_file
# Check if the earthquake is cleady in the case and warn the user if so
fle
=
h5py
.
File
(
self
.
cache_file
,
"a"
)
if
earthquake
.
id
in
list
(
fle
):
warnings
.
warn
(
"Earthquake %s already in cache file %s - "
"Running the shakemaps will overwrite this"
%
(
self
.
earthquake
.
id
,
self
.
cache_file
)
)
fle
.
close
()
else
:
self
.
caching
=
False
self
.
cache_file
=
None
self
.
rctx
=
None
self
.
dctx
=
None
self
.
sctx
=
None
self
.
synth_rupture_max_site_dist
=
valid
.
positive_float
(
synthetic_rupture_max_site_distance
,
"Max. Synthetic Rupture Site Distance"
)
self
.
synth_rupture_site_spacing
=
valid
.
positive_float
(
synthetic_rupture_site_spacing
,
"Synthetic Rupture Site Spacing"
)
self
.
_build_contexts
()
def
_build_contexts
(
self
):
"""
Construct the rupture, site and distances contexts from the earthquake, site and
ground motion models
"""
cmaker
=
ContextMaker
(
self
.
tectonic_region
,
self
.
ground_motion_model
[
"GMMs"
])
if
not
self
.
earthquake
.
rupture
:
# Use the finite rupture sampler to generate the rupture and corresponding
# distances from the available information
self
.
earthquake
.
_rupture
=
FiniteRuptureSampler
().
get_finite_rupture
(
self
.
num_rupture_samples
,
self
.
earthquake
,
rdim
=
self
.
rdim
,
weights
=
self
.
synthetic_rupture_distance_weights
,
maximum_site_distance
=
self
.
synth_rupture_max_site_dist
,
site_spacing
=
self
.
synth_rupture_site_spacing
,
)[
0
]
# For the sites and rupture context we can use the context maker to get all of the
# source, site and rupture distances
self
.
rctx
,
self
.
sctx
,
self
.
dctx
=
cmaker
.
make_contexts
(
self
.
site_model
.
get_site_collection
(),
self
.
earthquake
.
rupture
)
def
_cache_contexts
(
self
,
grp
:
h5py
.
Group
):
"""
If caching the shakemaps, then this stores the context information to the file
"""
# Add the contexts to the group object
ctxt
=
grp
.
create_group
(
"contexts"
)
rup_ctxt
=
ctxt
.
create_group
(
"rupture"
)
dist_ctxt
=
ctxt
.
create_group
(
"distances"
)
for
gmm
in
self
.
ground_motion_model
[
"GMMs"
]:
for
rup_attr
in
gmm
.
REQUIRES_RUPTURE_PARAMETERS
:
if
rup_attr
not
in
rup_ctxt
.
attrs
:
rup_ctxt
.
attrs
[
rup_attr
]
=
getattr
(
self
.
rctx
,
rup_attr
)
for
attr
in
gmm
.
REQUIRES_DISTANCES
:
if
attr
not
in
list
(
dist_ctxt
):
distance
=
getattr
(
self
.
dctx
,
attr
)
dist_dset
=
dist_ctxt
.
create_dataset
(
attr
,
distance
.
shape
,
dtype
=
"f"
)
dist_dset
[:]
=
distance
site_ctxt
=
ctxt
.
create_dataset
(
"sites"
,
self
.
site_model
.
site_array
.
shape
,
dtype
=
self
.
site_model
.
site_array
.
dtype
)
site_ctxt
[:]
=
self
.
site_model
.
site_array
if
self
.
site_model
.
bbox_properties
:
# If the site model has bounding box properties then store these
# as attributes
site_ctxt
.
attrs
[
"has_bbox"
]
=
True
site_ctxt
.
attrs
[
"llon"
]
=
self
.
site_model
.
bbox_properties
[
"bbox"
][
0
]
site_ctxt
.
attrs
[
"llat"
]
=
self
.
site_model
.
bbox_properties
[
"bbox"
][
1
]
site_ctxt
.
attrs
[
"ulon"
]
=
self
.
site_model
.
bbox_properties
[
"bbox"
][
2
]
site_ctxt
.
attrs
[
"ulat"
]
=
self
.
site_model
.
bbox_properties
[
"bbox"
][
3
]
site_ctxt
.
attrs
[
"spcx"
]
=
self
.
site_model
.
bbox_properties
[
"spcx"
]
site_ctxt
.
attrs
[
"spcy"
]
=
self
.
site_model
.
bbox_properties
[
"spcy"
]
site_ctxt
.
attrs
[
"ncol"
]
=
self
.
site_model
.
bbox_properties
[
"ncol"
]
site_ctxt
.
attrs
[
"nrow"
]
=
self
.
site_model
.
bbox_properties
[
"nrow"
]
else
:
site_ctxt
.
attrs
[
"has_bbox"
]
=
False
return
def
get_shakemap
(
self
,
intensity_measure_types
:
List
)
->
Tuple
[
np
.
ndarray
,
np
.
ndarray
,
Dict
]:
"""
Main function to constructs the shakemaps for the specified intensity measures,
caching the information to hdf5 if requested
Args:
intensity_measure_types: List of intensity measures for which the shakemaps will
be calculated
Returns:
aggregated_means: The mean of the ground motions from the different ground motion
models weighted by the assigned input weights
aggregated_stddevs: The total standard deviation of the ground motions from the
different ground motion models weighted by the assigned
input weights
shakemaps: Dictionary of individual shakemaps for each ground motion model
"""
shakemaps
=
{}
shakemap_dtypes
=
np
.
dtype
(
[
(
intensity_measure_type
,
np
.
float64
)
for
intensity_measure_type
in
intensity_measure_types
]
)
if
self
.
caching
:
# If caching, open (or create) the file and cache the contexts
fle
=
h5py
.
File
(
self
.
cache_file
,
"r+"
)
if
self
.
earthquake
.
id
in
list
(
fle
):
del
fle
[
self
.
earthquake
.
id
]
fle_eq
=
fle
.
create_group
(
self
.
earthquake
.
id
)
self
.
_cache_contexts
(
fle_eq
)
# Pre-allocate the aggregated shakemaps with zeros
aggregated_means
=
np
.
zeros
(
self
.
site_model
.
shape
,
dtype
=
shakemap_dtypes
)
aggregated_stddevs
=
np
.
zeros
(
self
.
site_model
.
shape
,
dtype
=
shakemap_dtypes
)
for
weight
,
gmm
in
zip
(
self
.
ground_motion_model
[
"weights"
],
self
.
ground_motion_model
[
"GMMs"
]
):
gmm_str
=
gmm
.
__class__
.
__name__
shakemaps
[
gmm_str
]
=
{
"weight"
:
weight
,
"mean"
:
np
.
zeros
(
self
.
site_model
.
shape
,
dtype
=
shakemap_dtypes
),
"stddev"
:
np
.
zeros
(
self
.
site_model
.
shape
,
dtype
=
shakemap_dtypes
),
}
for
intensity_measure_type
in
intensity_measure_types
:
input_imt
=
imt
.
from_string
(
intensity_measure_type
)
try
:
mean
,
[
stddev
]
=
gmm
.
get_mean_and_stddevs
(
self
.
sctx
,
self
.
rctx
,
self
.
dctx
,
input_imt
,
[
const
.
StdDev
.
TOTAL
]
)
except
KeyError
:
warnings
.
warn
(
"Ground motion model %s not defined for intensity "
"measure type %s"
%
(
str
(
gmm
),
intensity_measure_type
)
)
continue
aggregated_means
[
intensity_measure_type
]
+=
weight
*
mean
aggregated_stddevs
[
intensity_measure_type
]
+=
weight
*
stddev
shakemaps
[
gmm_str
][
"mean"
][
intensity_measure_type
]
=
mean
shakemaps
[
gmm_str
][
"stddev"
][
intensity_measure_type
]
=
stddev
if
self
.
caching
:
self
.
_cache_shakemap
(
fle_eq
,
shakemaps
,
aggregated_means
,
aggregated_stddevs
,
shakemap_dtypes
)
fle
.
close
()
return
aggregated_means
,
aggregated_stddevs
,
shakemaps
def
_cache_shakemap
(
self
,
fle
:
h5py
.
Group
,
shakemaps
:
Dict
,
aggregated_means
:
np
.
ndarray
,
aggregated_stddevs
:
np
.
ndarray
,
shakemap_dtypes
:
np
.
dtype
,
):
"""
If caching is required then the shakemaps are written to the hdf5 file
Args:
fle: HDF5 group object to store the shakemaps
shakemaps: Dictionary of individual shakemaps for each ground motion model
aggregated_means: The mean of the ground motions from the different ground motion
models weighted by the assigned input weights
aggregated_stddevs: The total standard deviation of the ground motions from the
different ground motion models weighted by the assigned
input weights
shakemap_dtypes: IMT-dependent data type of the shakemaps
"""
shakemap_grp
=
fle
.
create_group
(
"shakemaps"
)
for
gmm_string
in
shakemaps
:
gmm_grp
=
shakemap_grp
.
create_group
(
gmm_string
)
gmm_grp
.
attrs
[
"weight"
]
=
shakemaps
[
gmm_string
][
"weight"
]
mean_dset
=
gmm_grp
.
create_dataset
(
"mean"
,
shakemaps
[
gmm_string
][
"mean"
].
shape
,
dtype
=
shakemap_dtypes
)
mean_dset
[:]
=
shakemaps
[
gmm_string
][
"mean"
]
stddev_dset
=
gmm_grp
.
create_dataset
(
"stddev"
,
shakemaps
[
gmm_string
][
"stddev"
].
shape
,
dtype
=
shakemap_dtypes
)
stddev_dset
[:]
=
shakemaps
[
gmm_string
][
"stddev"
]
# Store the agregated results
aggregated_grp
=
fle
.
create_group
(
"aggregated"
)
aggregated_mean_dset
=
aggregated_grp
.
create_dataset
(
"mean"
,
aggregated_means
.
shape
,
shakemap_dtypes
)
aggregated_mean_dset
[:]
=
aggregated_means
aggregated_stddev_dset
=
aggregated_grp
.
create_dataset
(
"stddev"
,
aggregated_stddevs
.
shape
,
shakemap_dtypes
)
aggregated_stddev_dset
[:]
=
aggregated_stddevs
shakyground2/site_model.py
View file @
cd06865f
...
...
@@ -5,7 +5,7 @@ from __future__ import annotations
import
warnings
import
numpy
as
np
import
pandas
as
pd
from
typing
import
Tuple
,
Optional
,
Dict
from
typing
import
Tuple
,
Optional
,
Dict
,
Union
from
openquake.hazardlib.site
import
SiteCollection
...
...
@@ -56,6 +56,14 @@ class SiteModel(object):
# Returns the number of sites in the site array
return
self
.
site_array
.
shape
[
0
]
def
__getitem__
(
self
,
key
:
Union
[
str
,
int
]):
# Returns the named column from the site array (if provided with a string),
# or the corresponding row if provided with an integer
if
isinstance
(
key
,
str
)
or
isinstance
(
key
,
int
):
return
self
.
site_array
[
key
]
else
:
raise
KeyError
(
"Site model has no attribute %s"
%
key
)
def
__repr__
(
self
):
# Returns an informative string regarding the composition of the sites,
# including the bounding box information when relevant
...
...
@@ -72,6 +80,10 @@ class SiteModel(object):
info_string
+=
bbox_string
return
info_string
@
property
def
shape
(
self
):
return
self
.
site_array
.
shape
[
0
]
@
property
def
dataframe
(
self
):
# Return the site array as a dataframe
...
...
shakyground2/synthetic_rupture_generator.py
0 → 100644
View file @
cd06865f
"""
Generates a set of synthetic ruptures based on a point source configuration
"""
import
numpy
as
np
from
typing
import
List
,
Tuple
,
Optional
from
scipy.stats
import
truncnorm
from
openquake.hazardlib.geo
import
Point
,
Mesh
from
openquake.hazardlib.geo.surface.base
import
BaseSurface
from
openquake.hazardlib.source.rupture
import
ParametricProbabilisticRupture
from
shakyground2.earthquake
import
Earthquake
from
shakyground2.site_model
import
SiteModel
# The hypocentre distributions are defined from the Next Generation Attenuation
# (NGA) West 2 rupture database
NGAWEST_HYPO_DISTRIBUTIONS
=
{
"All"
:
{
"Mean"
:
np
.
array
([
0.46851357
,
0.67460261
]),
"Sigma"
:
np
.
array
([
0.20691535
,
0.22040932
]),
"COV"
:
np
.
array
([[
0.04293594
,
0.00103987
],
[
0.00103987
,
0.04871868
]]),
"DipRange"
:
[
30.0
,
90.0
],
},
"SS"
:
{
"Mean"
:
np
.
array
([
0.48545206
,
0.64942746
]),
"Sigma"
:
np
.
array
([
0.22415657
,
0.21677068
]),
"COV"
:
np
.
array
([[
0.05064814
,
-
0.00603003
],
[
-
0.00603003
,
0.04736544
]]),
"DipRange"
:
[
80.0
,
90.0
],
},
"R"
:
{
"Mean"
:
np
.
array
([
0.4674859
,
0.58483914
]),
"Sigma"
:
np
.
array
([
0.16275562
,
0.22017015
]),
"COV"
:
np
.
array
([[
0.02673021
,
0.0113362
],
[
0.0113362
,
0.04891558
]]),
"DipRange"
:
[
25.0
,
50.0
],
},
"N"
:
{
"Mean"
:
np
.
array
([
0.50887254
,
0.82404
]),
"Sigma"
:
np
.
array
([
0.22416128
,
0.13647917
]),
"COV"
:
np
.
array
([[
0.05085368
,
-
0.00332741
],
[
-
0.00332741
,
0.01885098
]]),
"DipRange"
:
[
50.0
,
80.0
],
},
}
# In case other definitions (thrust fault, undefined etc.) are used
NGAWEST_HYPO_DISTRIBUTIONS
[
"TF"
]
=
NGAWEST_HYPO_DISTRIBUTIONS
[
"R"
]
NGAWEST_HYPO_DISTRIBUTIONS
[
"U"
]
=
NGAWEST_HYPO_DISTRIBUTIONS
[
"All"
]