# ------------------- The OpenQuake Model Building Toolkit --------------------
# ------------------- FERMI: Fault nEtwoRks ModellIng -------------------------
# Copyright (C) 2023 GEM Foundation
# .-.
# / \ .-.
# | .`. ; .--. ___ .-. ___ .-. .-. ( __)
# | |(___) / \ ( ) \ ( ) ' \ (''")
# | |_ | .-. ; | ' .-. ; | .-. .-. ; | |
# ( __) | | | | | / (___) | | | | | | | |
# | | | |/ | | | | | | | | | | |
# | | | ' _.' | | | | | | | | | |
# | | | .'.-. | | | | | | | | | |
# | | ' `-' / | | | | | | | | | |
# (___) `.__.' (___) (___)(___)(___)(___)
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Affero General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more
# details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# -----------------------------------------------------------------------------
# vim: tabstop=4 shiftwidth=4 softtabstop=4
# coding: utf-8
import numpy as np
import numpy.typing as npt
from numba import njit, int64, boolean
from openquake.fnm.msr import area_to_mag
from openquake.hazardlib.geo.surface import SimpleFaultSurface
[docs]
def get_mags_and_areas(rups: list, areas: npt.ArrayLike, key: str = "generic"):
"""
Computes for each rupture the corresponding magnitude
:param rups:
A list of list with the indexes of the single-section ruptures forming
each rupture
:param areas:
An array with the area for each single section rupture.
:param key:
A key identifying the magnitude scaling relationship to be used.
:returns:
A :class:`numpy.ndarray` instance with the same cardinality of `rups`
"""
mags = np.zeros(len(rups))
tot_areas = np.zeros(len(rups))
for i_rup, rup in enumerate(rups):
totarea = np.sum([areas[i] for i in rup])
mags[i_rup] = area_to_mag(totarea, key)
tot_areas[i_rup] = totarea
return mags, tot_areas
def _get_rupture_area(surfs: list, rups: npt.ArrayLike) -> npt.ArrayLike:
"""
Computes the area for a single multi-fault rupture
:param surfs:
A list of surfaces i.e.
:class:`openquake.hazardlib.geo.surface.KiteSurfaces`
:param rups:
A :class:`numpy.ndarray` instance with the description of the rupture
:returns:
A float defining the area of the rupture
"""
area = 0.0
for subr in rups:
# Get the surface of the section
surf = surfs[int(subr[6])]
# Compute the surface of each cell composing the surface
#if hasattr(surf, "get_area"):
if isinstance(surf, SimpleFaultSurface):
surf_area = surf.get_area()
area += surf_area
else:
_, _, _, cell_area = surf.get_cell_dimensions()
idx = np.isfinite(cell_area)
i_row = np.arange(0, surf.mesh.lons.shape[0] - 1, 1)
i_col = np.arange(0, surf.mesh.lons.shape[1] - 1, 1)
mesh_col, mesh_row = np.meshgrid(i_col, i_row)
mask_row = np.logical_and(
mesh_row >= subr[0], mesh_row < subr[0] + subr[3]
)
mask_col = np.logical_and(
mesh_col >= subr[1], mesh_col < subr[1] + subr[2]
)
mask_all = np.logical_and(mask_row, mask_col)
mask = np.logical_and(mask_all, idx)
area += np.sum(cell_area[mask])
return area
[docs]
def get_ruptures_area(
surfs: npt.ArrayLike, rups: npt.ArrayLike
) -> npt.ArrayLike:
"""
Computes the area of single-section ruptures
:params surfs:
A :class:`numpy.ndarray` instance with the surfaces describing the
geometry of each section.
:params rups:
A :class:`numpy.ndarray` instance with the description of all the
single section ruptures.
:returns:
A :class:`numpy.ndarray` instance with the area [km2] of each single
single section rupture.
"""
areas = []
for rup in rups:
areas.append(_get_rupture_area(surfs, [rup]))
return np.array(areas)
def _get_ruptures_first_level(fault_system, aratios) -> npt.ArrayLike:
# Returns a :class:`numpy.ndarray` instance with the list of ruptures
# generated by each section in the fault system
#
# params:
# - fault_system
# - aratios
rups = []
nrups = 0
for i, row in enumerate(fault_system):
tmp = get_ruptures_section(
row[1], aspect_ratios=aratios, from_idx_rupture=nrups, idx=i
)
rups.extend(tmp)
nrups = len(rups)
rups = np.array(rups)
return rups
# @jit(int64[:](int64[:, :], int64[:], boolean[:], boolean[:]))
def _check_one_subrupture_has_connection(conns, subrupt, in_1st, in_2nd):
# Given:
# - `conns` the array with the connections between sections,
# - `subrupt` a rupture (or subrupture)
# - `in_1st` a boolean vector indicating if the connection for the current
# subrupture is the first one
# - `in_2nd` a boolean vector indicating if the connection for the current
# subrupture is the second one
#
# This function returns a tuple with the first element containing:
# - A boolean indicating if the subrupture contains the connection
# - An index that's 1 when the connection
# Define the range in columns and rows for the subrupture (in cells)
range_rows = np.array([subrupt[0], subrupt[0] + subrupt[3]], dtype=int)
range_cols = np.array([subrupt[1], subrupt[1] + subrupt[2]], dtype=int)
is_in = []
counter = 0
for i_conn in np.arange(0, conns.shape[0]):
conn = conns[i_conn, :]
if in_1st[i_conn]:
chk_row = np.logical_and(
conn[2] >= range_rows[0], conn[2] + conn[5] <= range_rows[1]
)
chk_col = np.logical_and(
conn[3] >= range_cols[0], conn[3] + conn[4] <= range_cols[1]
)
elif in_2nd[i_conn]:
chk_row = np.logical_and(
conn[6] >= range_rows[0], conn[6] + conn[9] <= range_rows[1]
)
chk_col = np.logical_and(
conn[7] >= range_cols[0], conn[7] + conn[8] <= range_cols[1]
)
else:
is_in.append([False, -1, -1, -1])
continue
conn_is_in = np.logical_and(chk_row, chk_col)
if not conn_is_in:
is_in.append([False, -1, -1, -1])
continue
# Index of the `other` section containing the subrupture
i_sec = conn[1] if in_1st[i_conn] else conn[0]
idx = 0 if in_1st[i_conn] else 1
# The last index will be updated in the calling function. `idx` is 1
# when the connection is the first one
is_in.append([conn_is_in, i_sec, idx, i_conn])
counter += 1
output = np.array(is_in, dtype=int)
return output
# @jit(int64[:](int64[:, :], int64[:]))
def _check_rupture_has_connections(connections, subrupt):
# Check the connections within a given rupture `subrupt`. Returns an array
# where each row contains the following information:
# - A boolean indicating if the subrupt contains a given connection
# - The index of the other section connected
# - A boolean. When True the other component of the connection is the
# first one provided (otherwise it's the second one)
# - The connection index (incremental). Can be used to select connections
# from the initial `connection` array
# Find the connections involving the section ID of the investigated
# sub-rupture
id_sect = int(subrupt[6])
cond_1 = connections[:, 0] == id_sect
cond_2 = connections[:, 1] == id_sect
idxs = np.where(np.logical_or(cond_1, cond_2))[0]
sub_conn = np.array(connections[idxs], dtype=int)
# Find if the subrupture in question contains any of the selected
# subsections
inp_1 = cond_1[idxs]
inp_2 = cond_2[idxs]
if sub_conn.ndim < 2:
sub_conn = np.expand_dims(sub_conn, axis=0)
tmp = _check_one_subrupture_has_connection(sub_conn, subrupt, inp_1, inp_2)
# Replacing the index for the subset of connections with the index of the
# connections found
if tmp.ndim > 1:
tmp[:, -1] = idxs[tmp[:, -1]]
return np.array(tmp, dtype=int)
# @jit(boolean(int64[:, :], int64[:, :]))
[docs]
def check_rup_exists(rups, srup):
"""
This function checks if `srup` is already included in the list of
ruptures `rups`
:param rups:
The set of ruptures
:param srup:
The rupture to be searched
"""
if rups.shape[1] < 1:
return False
# Check the indexes of subruptures with the same section and level
# and upper left corner
subr = srup[0, :]
idx_subrup = np.where(
(
(subr[0] == rups[:, 0])
& (subr[1] == rups[:, 1])
& (subr[6] == rups[:, 6])
& (subr[5] == rups[:, 5])
)
)
# If we find at least one match
check = False
if np.size(idx_subrup) > 0:
# Find indexes of ruptures at the same level (i.e. involving the same
# set of sections)
idx_rup = np.unique(rups[idx_subrup, 4])
# Create a subset of ruptures
tmp = np.where(np.isin(rups[:, 4], idx_rup))
rups_set = rups[tmp, :]
# Adjust the size of the array with the rupture
if rups_set.ndim > 2:
rups_set = np.squeeze(rups_set)
# Check if the subset contains the rupture
check = check_rup_exists_detail(rups_set, srup)
return check
[docs]
@njit(boolean(int64[:, :], int64[:, :]))
def check_rup_exists_detail(rups, srup):
"""
This function checks if `srup` is already included in the list of
ruptures `rups`
:param rups:
The set of ruptures
:param srup:
The rupture to be searched
"""
if rups.shape[1] < 1:
return False
# Process all the ruptures
for i_rup in np.unique(rups[:, 4]):
# List where we store the results of the match
match_subr = []
# Loop though the subruptures forming rupture with index `i_rup`
r_sub = rups[rups[:, 4] == i_rup, :]
for i_subr in np.arange(0, len(r_sub)):
# Get the subruptures
subr = r_sub[i_subr]
sr = srup[i_subr, :]
mask = np.array([0, 1, 2, 3, 5, 6])
chk = np.array_equal(subr[mask], sr[mask])
# If at least one is true it means that this subrupture is in
# the rupture
if chk:
match_subr.append(True)
else:
match_subr.append(False)
break
if np.all(np.array(match_subr)):
return True
return False
def _add_rups(srups, rup_sset, tmp, clevel, new_rup_idx, new_rups, i_rup):
# Create new ruptures by adding to the reference rupture
# the ones connected on the other section
for srup in srups[tmp]:
# Adding to the new rupture the subruptures of the
# reference rupture
new_rup = []
for subr in rup_sset[rup_sset[:, 4] == i_rup, :]:
subr[4] = new_rup_idx
subr[5] = clevel
new_rup.append(subr)
srup[4] = new_rup_idx
srup[5] = clevel
# Adding the current subrupture
new_rup.append(srup)
new_rup = np.array(new_rup)
# Fixing dimensions
if new_rup.ndim < 2:
new_rup = np.expand_dims(new_rup, axis=0)
tmp_rups = np.array(new_rups, dtype=int)
if tmp_rups.ndim < 2:
tmp_rups = np.expand_dims(tmp_rups, axis=0)
# Check
exists = check_rup_exists(tmp_rups, new_rup)
# Add the new rupture
if not exists:
new_rups.extend(new_rup)
new_rup_idx += 1
# @jit
[docs]
def get_ruptures_section(
ul_idx: np.ndarray,
aspect_ratios: np.ndarray = np.array([]),
from_idx_rupture: int = 0,
idx: int = -1,
) -> np.ndarray:
"""
Computes the ruptures admitted by a single section represented by:
the upper-left corner of each subsection and the subsection geometry
i.e. number of cells along strike and dip.
:param ul_idx:
An array with the upper left corner of each subsection
:param nc_stk_dip:
A vector with the lenght and width of subsections (in terms of
number of cells)
:param aspect_ratios:
An array with the minimum (included) and maximum value admitted
:param idx:
The index of the section generating these ruptures
:returns:
A :class:`np.ndarray` with the ruptures. Columns:
- 0: row of UL vertex
- 1: column of the UL vertex
- 2: number of cells (columns) forming the rupture
- 3: number of cells (rows) forming the rupture
- 4: index of the rupture
- 5: number of connections i.e. the level
- 6: the section ID
- 7: index of the rupture in this section
"""
# Set aspect ratios
if aspect_ratios.size < 1:
aspect_ratios = np.array([0.0, 1.0e10])
# This defines a rough estimate of the output
tmp = ul_idx.shape[0] ** 2 * ul_idx.shape[1] ** 2
# This defines the number of columns in the output matrix
shape_1 = 8
ruptures = np.zeros((tmp, shape_1), dtype=np.float64)
# Compute ruptures
c = 0
for irow in np.arange(0, ul_idx.shape[0]):
for icol in np.arange(0, ul_idx.shape[1]):
tmpa = [ul_idx[i, icol, 3] for i in np.arange(0, irow)]
tmpb = [ul_idx[irow, i, 2] for i in np.arange(0, icol)]
ruptures[c, 1] = np.sum(tmpb)
ruptures[c, 0] = np.sum(tmpa)
ruptures[c, 2] = ul_idx[irow, icol, 2]
ruptures[c, 3] = ul_idx[irow, icol, 3]
# Rupture incremental index for the current section
ruptures[c, 7] = c
c += 1
for rup_c in np.arange(icol, ul_idx.shape[1]):
for rup_r in np.arange(irow, ul_idx.shape[0]):
if rup_r == irow and rup_c == icol:
continue
ruptures[c, 0] = np.sum(tmpa)
ruptures[c, 1] = np.sum(tmpb)
tmp = [
ul_idx[i, icol, 3]
for i in np.arange(irow + 1, rup_r + 1)
]
ruptures[c, 3] = ul_idx[irow, icol, 3] + np.sum(tmp)
tmp = [
ul_idx[irow, i, 2]
for i in np.arange(icol + 1, rup_c + 1)
]
ruptures[c, 2] = ul_idx[irow, icol, 2] + np.sum(tmp)
ruptures[c, 7] = c
c += 1
# Remove empty rows
rups = np.array(ruptures[:c, :], dtype=np.float64)
# Set the level of the rupture
rups[:, 5] = 1.0
# Set the index of the section
if idx >= 0:
rups[:, 6] = idx
# Filter ruptures by their aspect ratio
asr = rups[:, 2] / rups[:, 3]
idx = (asr >= aspect_ratios[0]) & (asr < aspect_ratios[1])
tmp = np.arange(0, np.sum((idx)))
rups[idx, 4] = tmp + from_idx_rupture
rups = rups[idx, :]
# Fix the numbering of ruptures
rups[:, 7] = np.arange(0, rups.shape[0])
return rups