Coverage for /usr/local/lib/python3.11/dist-packages/grond/clustering/metrics.py: 32%
103 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-25 08:34 +0000
1import math
2import numpy as num
3from pyrocko import moment_tensor, orthodrome
4from grond.meta import GrondError
7def get_distance_mt_l2(eventi, eventj):
8 '''
9 L2 norm among two moment tensors, with 6 independet entries
10 '''
12 # ni = (eventi.mxx)**2 + (eventi.myy)**2 + (eventi.mzz)**2 + \
13 # (eventi.mxy)**2 + (eventi.mxz)**2 + (eventi.myz)**2
14 # nj = (eventj.mxx)**2 + (eventj.myy)**2 + (eventj.mzz)**2 + \
15 # (eventj.mxy)**2 + (eventj.mxz)**2 + (eventj.myz)**2
17 # mixx, miyy, mizz = eventi.mxx / ni, eventi.myy / ni, eventi.mzz / ni
18 # mixy, mixz, miyz = eventi.mxy / ni, eventi.mxz / ni, eventi.myz / ni
19 # mjxx, mjyy, mjzz = eventj.mxx / nj, eventj.myy / nj, eventj.mzz / nj
20 # mjxy, mjxz, mjyz = eventj.mxy / nj, eventj.mxz / nj, eventj.myz / nj
22 d = (eventi.mxx - eventj.mxx)**2 + (eventi.myy - eventj.myy)**2 + \
23 (eventi.mzz - eventj.mzz)**2 + (eventi.mxy - eventj.mxy)**2 + \
24 (eventi.mxz - eventj.mxz)**2 + (eventi.myz - eventj.myz)**2
26 d = 0.5 * math.sqrt(d)
28 return d
31def get_distance_mt_l1(eventi, eventj):
32 '''
33 L1 norm among two moment tensors, with 6 independet entries
34 '''
36 # ni = abs(eventi.mxx) + abs(eventi.myy) + abs(eventi.mzz) + \
37 # abs(eventi.mxy) + abs(eventi.mxz) + abs(eventi.myz)
38 # nj = abs(eventj.mxx) + abs(eventj.myy) + abs(eventj.mzz) + \
39 # abs(eventj.mxy) + abs(eventj.mxz) + abs(eventj.myz)
40 #
41 # mixx, miyy, mizz = eventi.mxx / ni, eventi.myy / ni, eventi.mzz / ni
42 # mixy, mixz, miyz = eventi.mxy / ni, eventi.mxz / ni, eventi.myz / ni
43 # mjxx, mjyy, mjzz = eventj.mxx / nj, eventj.myy / nj, eventj.mzz / nj
44 # mjxy, mjxz, mjyz = eventj.mxy / nj, eventj.mxz / nj, eventj.myz / nj
46 d = abs(eventi.mxx - eventj.mxx) + abs(eventi.myy - eventj.myy) + \
47 abs(eventi.mzz - eventj.mzz) + abs(eventi.mxy - eventj.mxy) + \
48 abs(eventi.mxz - eventj.mxz) + abs(eventi.myz - eventj.myz)
50 d = 0.5 * math.sqrt(d)
52 return d
55def get_distance_mt_cos(eventi, eventj):
56 '''
57 Inner product among two moment tensors.
59 According to Willemann 1993; and later to Tape & Tape, normalization in
60 R^9 to ensure innerproduct between -1 and +1.
61 '''
63 ni = math.sqrt(
64 eventi.mxx * eventi.mxx +
65 eventi.myy * eventi.myy +
66 eventi.mzz * eventi.mzz +
67 2. * eventi.mxy * eventi.mxy +
68 2. * eventi.mxz * eventi.mxz +
69 2. * eventi.myz * eventi.myz)
71 nj = math.sqrt(
72 eventj.mxx * eventj.mxx +
73 eventj.myy * eventj.myy +
74 eventj.mzz * eventj.mzz +
75 2. * eventj.mxy * eventj.mxy +
76 2. * eventj.mxz * eventj.mxz +
77 2. * eventj.myz * eventj.myz)
79 nc = ni * nj
80 innerproduct = (
81 eventi.mxx * eventj.mxx +
82 eventi.myy * eventj.myy +
83 eventi.mzz * eventj.mzz +
84 2. * eventi.mxy * eventj.mxy +
85 2. * eventi.mxz * eventj.mxz +
86 2. * eventi.myz * eventj.myz) / nc
88 if innerproduct >= 1.0:
89 innerproduct = 1.0
90 elif innerproduct <= -1.0:
91 innerproduct = -1.0
93 d = 0.5 * (1 - innerproduct)
95 return d
98def get_distance_mt_weighted_cos(eventi, eventj, ws):
99 '''
100 Weighted moment tensor distance.
102 According to Cesca et al. 2014 GJI
103 '''
104 ni = math.sqrt(
105 (ws[0] * eventi.mxx)**2 +
106 (ws[1] * eventi.mxy)**2 +
107 (ws[2] * eventi.myy)**2 +
108 (ws[3] * eventi.mxz)**2 +
109 (ws[4] * eventi.myz)**2 +
110 (ws[5] * eventi.mzz)**2)
112 nj = math.sqrt(
113 (ws[0] * eventj.mxx)**2 +
114 (ws[1] * eventj.mxy)**2 +
115 (ws[2] * eventj.myy)**2 +
116 (ws[3] * eventj.mxz)**2 +
117 (ws[4] * eventj.myz)**2 +
118 (ws[5] * eventj.mzz)**2)
120 nc = ni * nj
121 innerproduct = (
122 ws[0] * ws[0] * eventi.mxx * eventj.mxx +
123 ws[1] * ws[1] * eventi.mxy * eventj.mxy +
124 ws[2] * ws[2] * eventi.myy * eventj.myy +
125 ws[3] * ws[3] * eventi.mxz * eventj.mxz +
126 ws[4] * ws[4] * eventi.myz * eventj.myz +
127 ws[5] * ws[5] * eventi.mzz * eventj.mzz) / nc
129 if innerproduct >= 1.0:
130 innerproduct = 1.0
132 elif innerproduct <= -1.0:
133 innerproduct = -1.0
135 d = 0.5 * (1.0 - innerproduct)
137 return d
140def get_distance_dc(eventi, eventj):
141 '''Normalized Kagan angle distance among DC components of moment tensors.
142 Based on Kagan, Y. Y., 1991, GJI
143 '''
145 mti = eventi.moment_tensor
146 mtj = eventj.moment_tensor
148 d = moment_tensor.kagan_angle(mti, mtj) / 120.
149 if d > 1.:
150 d = 1.
152 return d
155def get_distance_hypo(eventi, eventj):
156 '''
157 Normalized Euclidean hypocentral distance, assuming flat earth to combine
158 epicentral distance and depth difference.
160 The normalization assumes largest considered distance is 1000 km.
161 '''
162 maxdist_km = 1000.
163 a_lats, a_lons, b_lats, b_lons = \
164 eventi.north, eventi.east, eventj.north, eventj.east
166 a_dep, b_dep = eventi.down, eventj.down
168 if (a_lats == b_lats) and (a_lons == b_lons) and (a_dep == b_dep):
169 d = 0.
170 else:
171 distance_m = orthodrome.distance_accurate50m_numpy(
172 a_lats, a_lons, b_lats, b_lons)
174 distance_km = distance_m / 1000.
175 ddepth = abs(eventi.down - eventj.down)
176 hypo_distance_km = math.sqrt(
177 distance_km * distance_km + ddepth * ddepth)
179 # maxdist = float(inv_param['EUCLIDEAN_MAX'])
181 d = hypo_distance_km / maxdist_km
182 if d >= 1.:
183 d = 1.
185 return d
188def get_distance_epi(eventi, eventj):
189 '''Normalized Euclidean epicentral distance.
190 The normalization assumes largest considered distance is 1000 km.
191 '''
192 maxdist_km = 1000.
194 a_lats, a_lons, b_lats, b_lons = \
195 eventi.north, eventi.east, eventj.north, eventj.east
197 a_dep, b_dep = eventi.down, eventj.down
199 if (a_lats == b_lats) and (a_lons == b_lons) and (a_dep == b_dep):
200 d = 0.
201 else:
202 distance_m = orthodrome.distance_accurate50m_numpy(
203 a_lats, a_lons, b_lats, b_lons)
205 distance_km = distance_m / 1000.
207 d = distance_km / maxdist_km
208 if d >= 1.:
209 d = 1.
211 return d
214def get_distance_mt_triangle_diagram(eventi, eventj):
215 '''
216 Scalar product among principal axes (?).
217 '''
219 mti = eventi.moment_tensor
220 mtj = eventj.moment_tensor
222 ti, pi, bi = mti.t_axis(), mti.p_axis(), mti.b_axis()
223 deltabi = math.acos(abs(bi[2]))
224 deltati = math.acos(abs(ti[2]))
225 deltapi = math.acos(abs(pi[2]))
226 tj, pj, bj = mtj.t_axis(), mtj.p_axis(), mtj.b_axis()
227 deltabj = math.acos(abs(bj[2]))
228 deltatj = math.acos(abs(tj[2]))
229 deltapj = math.acos(abs(pj[2]))
230 dotprod = deltabi * deltabj + deltati * deltatj + deltapi * deltapj
232 if dotprod >= 1.:
233 dotprod == 1.
235 d = 1. - dotprod
236 return d
239metric_funcs = {
240 'mt_l2norm': get_distance_mt_l2,
241 'mt_l1norm': get_distance_mt_l1,
242 'mt_cos': get_distance_mt_cos,
243 'mt_weighted_cos': get_distance_mt_weighted_cos,
244 'mt_principal_axis': get_distance_mt_triangle_diagram,
245 'kagan_angle': get_distance_dc,
246 'hypocentral': get_distance_hypo,
247 'epicentral': get_distance_epi,
248}
251metrics = sorted(metric_funcs.keys())
254def get_distance(eventi, eventj, metric, **kwargs):
255 '''
256 Compute the normalized distance among two earthquakes, calling the function
257 for the chosen metric definition.
258 '''
260 try:
261 func = metric_funcs[metric]
262 except KeyError:
263 raise GrondError('unknown metric: %s' % metric)
265 return func(eventi, eventj)
268def compute_similarity_matrix(events, metric):
269 '''
270 Compute and return a similarity matrix for all event pairs, according to
271 the desired metric
273 :param events: list of pyrocko events
274 :param metric: metric type (string)
276 :returns: similarity matrix as NumPy array
277 '''
279 nev = len(events)
280 simmat = num.zeros((nev, nev), dtype=float)
281 for i in range(len(events)):
282 for j in range(i):
283 d = get_distance(events[i], events[j], metric)
284 simmat[i, j] = d
285 simmat[j, i] = d
287 return simmat
290def load_similarity_matrix(fname):
291 '''
292 Load a binary similarity matrix from file
293 '''
295 simmat = num.load(fname)
296 return simmat