# ------------------- The OpenQuake Model Building Toolkit --------------------
# Copyright (C) 2022 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 os
import h5py
import numpy as np
import matplotlib.pyplot as plt
import pickle
import logging
from decimal import Decimal, getcontext
from scipy.interpolate import RBFInterpolator
from openquake.mbt.tools.tr.catalogue import get_catalogue
from openquake.mbt.tools.geo import get_idx_points_inside_polygon
from openquake.mbt.tools.tr.catalogue_hmtk import (get_rtree_index,
get_distances_from_surface)
from openquake.sub.utils import (_read_edges,
build_complex_surface_from_edges,
build_kite_surface_from_profiles,
plot_complex_surface)
from openquake.hmtk.seismicity.selector import CatalogueSelector
getcontext().prec = 4
# Buffer around the brounding box
#DELTA = 0.3
DELTA = 1
[docs]
class SetSubductionEarthquakes:
"""
Classifies earthquakes generated by a subduction earthquake source (either
interface or inslab).
:param str label:
Flag used to classify the earthquakes
:param str treg_filename:
Name of the .hdf5 containing the tectonic regionalisation
:param str distance_folder:
Folder where to store the epicenter-surface files
:param str edges_folder:
Folder containing the edges specifying the geometry of the surface
:param float distance_buffer_below:
Distance [km] below the fault surface used to select earthquakes
:param float distance_buffer_above:
Distance [km] above the fault surface used to select earthquakes
:param str catalogue_filename:
Name of the file containing the earthquakes to be classified
:param str log_fname:
Name of the .log file
:param low_year:
Lowest year selected
:param upp_year:
Largest year selected
:param low_mag:
Lowest magnitude selected
:param upp_mag:
Largest magnitude selected
"""
def __init__(self, label, treg_filename, distance_folder, edges_folder,
distance_buffer_below, distance_buffer_above, lower_depth,
catalogue_filename, log_fname=None, upper_depth=None,
low_year=-10000, upp_year=+10000, low_mag=-5., upp_mag=15.):
self.label = label
self.treg_filename = treg_filename
self.distance_folder = distance_folder
self.edges_folder = edges_folder
self.distance_buffer_below = Decimal(distance_buffer_below)
self.distance_buffer_above = Decimal(distance_buffer_above)
self.catalogue_filename = catalogue_filename
self.lower_depth = lower_depth
self.upper_depth = upper_depth
self.log_fname = log_fname
self.low_year = low_year
self.upp_year = upp_year
self.low_mag = Decimal(low_mag)
self.upp_mag = Decimal(upp_mag)
[docs]
def classify(self, compute_distances, remove_from, surftype='ComplexFault'):
"""
:param bool compute_distances:
A boolean indicating if distances between earthquakes and the
subduction surface should be computed. If False the distances
stored in `self.distance_folder` will be used.
:param list remove_from:
A list of labels identifying TR from where the earthquakes assigned
to this TR must be removed
"""
# set parameters
treg_filename = self.treg_filename
distance_folder = self.distance_folder
edges_folder = self.edges_folder
distance_buffer_below = self.distance_buffer_below
distance_buffer_above = self.distance_buffer_above
catalogue_filename = self.catalogue_filename
lower_depth = self.lower_depth
if lower_depth is None:
lower_depth = 400.
upper_depth = self.upper_depth
if upper_depth is None:
upper_depth = 0
# Open log file and prepare the group
print(f'Log filename: {self.log_fname}\n')
flog = h5py.File(self.log_fname, 'a')
if self.label not in flog.keys():
grp = flog.create_group('/{:s}'.format(self.label))
else:
grp = flog['/{:s}'.format(self.label)]
# Read the catalogue
catalogue = get_catalogue(catalogue_filename)
neq = len(catalogue.data['longitude'])
f = h5py.File(treg_filename, "a")
if self.label in f.keys():
treg = f[self.label]
else:
treg = np.full((neq), False, dtype=bool)
# Create the spatial index
sidx = get_rtree_index(catalogue)
# Build the complex fault surface
tedges = _read_edges(edges_folder)
print(edges_folder)
if surftype == 'ComplexFault':
surface = build_complex_surface_from_edges(edges_folder)
# Create polygon encompassing the mesh
mesh = surface.mesh
plo = list(mesh.lons[0, :])
pla = list(mesh.lats[0, :])
#
plo += list(mesh.lons[:, -1])
pla += list(mesh.lats[:, -1])
#
plo += list(mesh.lons[-1, ::-1])
pla += list(mesh.lats[-1, ::-1])
#
plo += list(mesh.lons[::-1, 0])
pla += list(mesh.lats[::-1, 0])
elif surftype == 'KiteFault':
surface = build_kite_surface_from_profiles(edges_folder)
# Create polygon encompassing the mesh
mesh = surface.mesh
plo = surface.surface_projection[0]
pla = surface.surface_projection[1]
else:
msg = f'surface type {surftype} not supported'
raise ValueError(msg)
# Set variables used in griddata
data = np.array([mesh.lons.flatten().T, mesh.lats.flatten().T]).T
values = mesh.depths.flatten().T
depths_dec = [Decimal(x) for x in values]
depths = np.array(depths_dec)
ddd = np.array([mesh.lons.flatten().T,
mesh.lats.flatten().T,
depths]).T
if self.label not in flog.keys():
grp.create_dataset('mesh', data=ddd)
# Set the bounding box of the subduction surface
min_lo_sub = np.nanmin(mesh.lons)
min_la_sub = np.nanmin(mesh.lats)
max_lo_sub = np.nanmax(mesh.lons)
max_la_sub = np.nanmax(mesh.lats)
# Select the earthquakes within the bounding box
idxs = sorted(list(sidx.intersection((min_lo_sub-DELTA,
min_la_sub-DELTA,
upper_depth,
max_lo_sub+DELTA,
max_la_sub+DELTA,
lower_depth))))
# Select earthquakes within the bounding box of the surface
# projection of the fault
sidx = get_idx_points_inside_polygon(catalogue.data['longitude'][idxs],
catalogue.data['latitude'][idxs],
plo, pla,
idxs, buff_distance=5000.)
# Select earthquakes and store indexes of selected ones
ccc = []
idxs = []
for idx in sidx:
# Preselection based on magnitude and time of occurrence
if ((catalogue.data['magnitude'][idx] >= self.low_mag) &
(catalogue.data['magnitude'][idx] <= self.upp_mag) &
(catalogue.data['year'][idx] >= self.low_year) &
(catalogue.data['year'][idx] <= self.upp_year)):
idxs.append(idx)
# Update the log file
ccc.append([catalogue.data['longitude'][idx],
catalogue.data['latitude'][idx],
catalogue.data['depth'][idx]])
if self.label not in flog.keys():
grp.create_dataset('cat', data=np.array(ccc))
# Prepare array for the selection of the catalogue
flags = np.full((len(catalogue.data['longitude'])), False, dtype=bool)
flags[idxs] = True
# Create a selector for the catalogue and select earthquakes within
# bounding box
sel = CatalogueSelector(catalogue, create_copy=True)
cat = sel.select_catalogue(flags)
self.cat = cat
# If none of the earthquakes in the catalogue is in the bounding box
# used for the selection we stop the processing
if len(cat.data['longitude']) < 1:
f = h5py.File(treg_filename, "a")
if self.label in f.keys():
del f[self.label]
f[self.label] = treg
f.close()
return
# compute distances between the earthquakes in the catalogue and
# the surface of the fault
out_filename = os.path.join(distance_folder,
'dist_{:s}.pkl'.format(self.label))
surf_dist = get_distances_from_surface(cat, surface)
"""
if compute_distances:
tmps = 'Computing distances'
logging.info(tmps.format(out_filename))
surf_dist = get_distances_from_surface(cat, surface)
pickle.dump(surf_dist, open(out_filename, 'wb'))
else:
if not os.path.exists(out_filename):
raise IOError('Distance file does not exist')
surf_dist = pickle.load(open(out_filename, 'rb'))
tmps = 'Loading distances from file: {:s}'
logging.info(tmps.format(out_filename))
tmps = ' number of values loaded: {:d}'
logging.info(tmps.format(len(surf_dist)))
"""
# info
neqks = len(cat.data['longitude'])
tmps = 'Number of eqks in the new catalogue : {:d}'
logging.info(tmps.format(neqks))
# Calculate the depth of the top of the slab for every earthquake
# location
points = np.array([[lo, la] for lo, la in zip(cat.data['longitude'],
cat.data['latitude'])])
# Compute the depth of the top of the slab at every epicenter using
# interpolation
# sub_depths = griddata(data, values, (points[:, 0], points[:, 1]),
# method='cubic')
val_red = values[~np.isnan(values)]
dat_red = data[~np.isnan(data)].reshape(-1, 2)
# dat_red_fi = dat_red.reshape(len(val_red), 2)
rbfi = RBFInterpolator(dat_red[:, 0:2], val_red, kernel='multiquadric',
epsilon=1, neighbors=100)
sub_depths = rbfi(points[:, 0:2])
# Save the distances to a file
tmps = 'vert_dist_to_slab_{:s}.pkl'.format(self.label)
out_filename = os.path.join(distance_folder, tmps)
if not os.path.exists(out_filename):
pickle.dump(surf_dist, open(out_filename, 'wb'))
# Let's find earthquakes close to the top of the slab
idxa = np.nonzero((np.isfinite(surf_dist) &
np.isfinite(sub_depths) &
np.isfinite(cat.data['depth'])) &
((surf_dist < distance_buffer_below) &
(sub_depths > cat.data['depth'])) |
((surf_dist < distance_buffer_above) &
(sub_depths <= cat.data['depth'])))[0]
idxa = []
for srfd, subd, dept in zip(surf_dist, sub_depths, cat.data['depth']):
if np.isfinite(srfd) & np.isfinite(subd) & np.isfinite(dept):
if (Decimal(srfd) < min(distance_buffer_below,
distance_buffer_above) * Decimal(0.90)):
idxa.append(True)
elif ((Decimal(srfd) < distance_buffer_below) &
(Decimal(subd) < Decimal(dept))):
idxa.append(True)
elif ((Decimal(srfd) < distance_buffer_above) &
(Decimal(subd) >= Decimal(dept))):
idxa.append(True)
else:
idxa.append(False)
else:
idxa.append(False)
idxa = np.array(idxa)
print(idxa)
# Check the size of lists
assert len(idxa) == len(cat.data['longitude']) == len(idxs)
self.surf_dist = surf_dist
self.sub_depths = sub_depths
self.tedges = tedges
self.idxa = idxa
self.treg = treg
# Store log data
tl = np.zeros(len(idxa),
dtype={'names': ('eid', 'lon', 'lat', 'dep', 'subd',
'srfd', 'idx'),
'formats': ('S15', 'f8', 'f8', 'f8', 'f8', 'f8',
'i4')})
tl['eid'] = cat.data['eventID']
tl['lon'] = cat.data['longitude']
tl['lat'] = cat.data['latitude']
tl['dep'] = cat.data['depth']
tl['subd'] = sub_depths
tl['srfd'] = surf_dist
tl['idx'] = idxa
if not 'data' in grp.keys():
grp.create_dataset('data', data=np.array(tl))
# Update the selection array
for uuu, iii in enumerate(list(idxa)):
aaa = idxs[uuu]
assert catalogue.data['eventID'][aaa] == cat.data['eventID'][uuu]
if iii:
treg[aaa] = True
else:
treg[aaa] = False
# Store results in the .hdf5 file
logging.info('Storing data in:\n{:s}'.format(treg_filename))
f = h5py.File(treg_filename, "a")
if len(remove_from):
fmt = ' treg: {:d}'
logging.info(fmt.format(len(treg)))
iii = np.nonzero(treg)[0]
for tkey in remove_from:
logging.info(' Cleaning {:s}'.format(tkey))
old = f[tkey][:]
fmt = ' before: {:d}'
logging.info(fmt.format(len(np.nonzero(old)[0])))
del f[tkey]
old[iii] = False
f[tkey] = old
fmt = ' after: {:d}'
logging.info(fmt.format(len(np.nonzero(old)[0])))
# Remove the old classification and adding the new one
if self.label in f.keys():
del f[self.label]
f[self.label] = treg
# Close files
f.close()
flog.close()
[docs]
def plotting_0(self):
"""
"""
cat = self.cat
sub_depths = self.sub_depths
surf_dist = self.surf_dist
#
plt.figure(figsize=(10, 8))
scat = plt.scatter(cat.data['depth'], sub_depths, c=surf_dist,
s=2**cat.data['magnitude'], edgecolor='w', vmin=0,
vmax=100)
plt.ylabel('Top of slab depth [km]')
plt.xlabel('Hypocentral depth [km]')
xx = np.arange(10, 300)
plt.plot(xx, xx, ':r')
plt.xscale('log')
plt.yscale('log')
plt.xlim([10, 300])
plt.ylim([10, 200])
plt.grid(axis='both', which='both')
cb = plt.colorbar(scat, extend='both')
cb.set_label('Shortest distance [km]')
[docs]
def plotting_1(self):
"""
"""
cat = self.cat
tedges = self.tedges
idxa = self.idxa
fig, ax = plot_complex_surface(tedges)
ax.plot(cat.data['longitude'],
cat.data['latitude'],
cat.data['depth'], '.b', alpha=0.05)
ax.plot(cat.data['longitude'][idxa],
cat.data['latitude'][idxa],
cat.data['depth'][idxa], '.c', alpha=0.2)
plt.show()