Coverage for /usr/local/lib/python3.11/dist-packages/grond/clustering/metrics.py: 31%

1# https://pyrocko.org/grond - GPLv3 

2# 

3# The Grond Developers, 21st Century 

4import math 

5import numpy as num 

6from pyrocko import moment_tensor, orthodrome 

7from grond.meta import GrondError 

8 

9 

10def get_distance_mt_l2(eventi, eventj): 

11 ''' 

12 L2 norm among two moment tensors, with 6 independet entries 

13 ''' 

14 

15 m1 = eventi.moment_tensor.m() 

16 m2 = eventj.moment_tensor.m() 

17 

18 m1 /= num.linalg.norm(m1) 

19 m2 /= num.linalg.norm(m2) 

20 

21 return 0.5 * num.sqrt(num.sum((m1 - m2)**2)) 

22 

23 

24def get_distance_mt_l1(eventi, eventj): 

25 ''' 

26 L1 norm among two moment tensors, with 6 independet entries 

27 ''' 

28 

29 m1 = eventi.moment_tensor.m() 

30 m2 = eventj.moment_tensor.m() 

31 

32 m1 /= num.linalg.norm(m1) 

33 m2 /= num.linalg.norm(m2) 

34 

35 return 0.5 * num.sum( 

36 num.abs(m1 - m2)) 

37 

38 

39def get_distance_mt_cos(eventi, eventj): 

40 ''' 

41 Inner product among two moment tensors. 

42 

43 According to Willemann 1993; and later to Tape & Tape, normalization in 

44 R^9 to ensure innerproduct between -1 and +1. 

45 ''' 

46 

47 ni = math.sqrt( 

48 eventi.mxx * eventi.mxx + 

49 eventi.myy * eventi.myy + 

50 eventi.mzz * eventi.mzz + 

51 2. * eventi.mxy * eventi.mxy + 

52 2. * eventi.mxz * eventi.mxz + 

53 2. * eventi.myz * eventi.myz) 

54 

55 nj = math.sqrt( 

56 eventj.mxx * eventj.mxx + 

57 eventj.myy * eventj.myy + 

58 eventj.mzz * eventj.mzz + 

59 2. * eventj.mxy * eventj.mxy + 

60 2. * eventj.mxz * eventj.mxz + 

61 2. * eventj.myz * eventj.myz) 

62 

63 nc = ni * nj 

64 innerproduct = ( 

65 eventi.mxx * eventj.mxx + 

66 eventi.myy * eventj.myy + 

67 eventi.mzz * eventj.mzz + 

68 2. * eventi.mxy * eventj.mxy + 

69 2. * eventi.mxz * eventj.mxz + 

70 2. * eventi.myz * eventj.myz) / nc 

71 

72 if innerproduct >= 1.0: 

73 innerproduct = 1.0 

74 elif innerproduct <= -1.0: 

75 innerproduct = -1.0 

76 

77 d = 0.5 * (1 - innerproduct) 

78 

79 return d 

80 

81 

82def get_distance_mt_weighted_cos(eventi, eventj, ws): 

83 ''' 

84 Weighted moment tensor distance. 

85 

86 According to Cesca et al. 2014 GJI 

87 ''' 

88 ni = math.sqrt( 

89 (ws[0] * eventi.mxx)**2 + 

90 (ws[1] * eventi.mxy)**2 + 

91 (ws[2] * eventi.myy)**2 + 

92 (ws[3] * eventi.mxz)**2 + 

93 (ws[4] * eventi.myz)**2 + 

94 (ws[5] * eventi.mzz)**2) 

95 

96 nj = math.sqrt( 

97 (ws[0] * eventj.mxx)**2 + 

98 (ws[1] * eventj.mxy)**2 + 

99 (ws[2] * eventj.myy)**2 + 

100 (ws[3] * eventj.mxz)**2 + 

101 (ws[4] * eventj.myz)**2 + 

102 (ws[5] * eventj.mzz)**2) 

103 

104 nc = ni * nj 

105 innerproduct = ( 

106 ws[0] * ws[0] * eventi.mxx * eventj.mxx + 

107 ws[1] * ws[1] * eventi.mxy * eventj.mxy + 

108 ws[2] * ws[2] * eventi.myy * eventj.myy + 

109 ws[3] * ws[3] * eventi.mxz * eventj.mxz + 

110 ws[4] * ws[4] * eventi.myz * eventj.myz + 

111 ws[5] * ws[5] * eventi.mzz * eventj.mzz) / nc 

112 

113 if innerproduct >= 1.0: 

114 innerproduct = 1.0 

115 

116 elif innerproduct <= -1.0: 

117 innerproduct = -1.0 

118 

119 d = 0.5 * (1.0 - innerproduct) 

120 

121 return d 

122 

123 

124def get_distance_dc(eventi, eventj): 

125 '''Normalized Kagan angle distance among DC components of moment tensors. 

126 Based on Kagan, Y. Y., 1991, GJI 

127 ''' 

128 

129 mti = eventi.moment_tensor 

130 mtj = eventj.moment_tensor 

131 

132 d = moment_tensor.kagan_angle(mti, mtj) / 120. 

133 if d > 1.: 

134 d = 1. 

135 

136 return d 

137 

138 

139def get_distance_hypo(eventi, eventj): 

140 ''' 

141 Normalized Euclidean hypocentral distance, assuming flat earth to combine 

142 epicentral distance and depth difference. 

143 

144 The normalization assumes largest considered distance is 1000 km. 

145 ''' 

146 maxdist_km = 1000. 

147 a_lats, a_lons, b_lats, b_lons = \ 

148 eventi.north, eventi.east, eventj.north, eventj.east 

149 

150 a_dep, b_dep = eventi.down, eventj.down 

151 

152 if (a_lats == b_lats) and (a_lons == b_lons) and (a_dep == b_dep): 

153 d = 0. 

154 else: 

155 distance_m = orthodrome.distance_accurate50m_numpy( 

156 a_lats, a_lons, b_lats, b_lons) 

157 

158 distance_km = distance_m / 1000. 

159 ddepth = abs(eventi.down - eventj.down) 

160 hypo_distance_km = math.sqrt( 

161 distance_km * distance_km + ddepth * ddepth) 

162 

163 # maxdist = float(inv_param['EUCLIDEAN_MAX']) 

164 

165 d = hypo_distance_km / maxdist_km 

166 if d >= 1.: 

167 d = 1. 

168 

169 return d 

170 

171 

172def get_distance_epi(eventi, eventj): 

173 '''Normalized Euclidean epicentral distance. 

174 The normalization assumes largest considered distance is 1000 km. 

175 ''' 

176 maxdist_km = 1000. 

177 

178 a_lats, a_lons, b_lats, b_lons = \ 

179 eventi.north, eventi.east, eventj.north, eventj.east 

180 

181 a_dep, b_dep = eventi.down, eventj.down 

182 

183 if (a_lats == b_lats) and (a_lons == b_lons) and (a_dep == b_dep): 

184 d = 0. 

185 else: 

186 distance_m = orthodrome.distance_accurate50m_numpy( 

187 a_lats, a_lons, b_lats, b_lons) 

188 

189 distance_km = distance_m / 1000. 

190 

191 d = distance_km / maxdist_km 

192 if d >= 1.: 

193 d = 1. 

194 

195 return d 

196 

197 

198def get_distance_mt_triangle_diagram(eventi, eventj): 

199 ''' 

200 Scalar product among principal axes (?). 

201 ''' 

202 

203 mti = eventi.moment_tensor 

204 mtj = eventj.moment_tensor 

205 

206 ti, pi, bi = mti.t_axis(), mti.p_axis(), mti.b_axis() 

207 deltabi = math.acos(abs(bi[2])) 

208 deltati = math.acos(abs(ti[2])) 

209 deltapi = math.acos(abs(pi[2])) 

210 tj, pj, bj = mtj.t_axis(), mtj.p_axis(), mtj.b_axis() 

211 deltabj = math.acos(abs(bj[2])) 

212 deltatj = math.acos(abs(tj[2])) 

213 deltapj = math.acos(abs(pj[2])) 

214 dotprod = deltabi * deltabj + deltati * deltatj + deltapi * deltapj 

215 

216 if dotprod >= 1.: 

217 dotprod == 1. 

218 

219 d = 1. - dotprod 

220 return d 

221 

222 

223metric_funcs = { 

224 'mt_l2norm': get_distance_mt_l2, 

225 'mt_l1norm': get_distance_mt_l1, 

226 'mt_cos': get_distance_mt_cos, 

227 'mt_weighted_cos': get_distance_mt_weighted_cos, 

228 'mt_principal_axis': get_distance_mt_triangle_diagram, 

229 'kagan_angle': get_distance_dc, 

230 'hypocentral': get_distance_hypo, 

231 'epicentral': get_distance_epi, 

232} 

233 

234 

235metrics = sorted(metric_funcs.keys()) 

236 

237 

238def get_distance(eventi, eventj, metric, **kwargs): 

239 ''' 

240 Compute the normalized distance among two earthquakes, calling the function 

241 for the chosen metric definition. 

242 ''' 

243 

244 try: 

245 func = metric_funcs[metric] 

246 except KeyError: 

247 raise GrondError('unknown metric: %s' % metric) 

248 

249 return func(eventi, eventj) 

250 

251 

252def compute_similarity_matrix(events, metric): 

253 ''' 

254 Compute and return a similarity matrix for all event pairs, according to 

255 the desired metric 

256 

257 :param events: list of pyrocko events 

258 :param metric: metric type (string) 

259 

260 :returns: similarity matrix as NumPy array 

261 ''' 

262 

263 nev = len(events) 

264 simmat = num.zeros((nev, nev), dtype=float) 

265 for i in range(len(events)): 

266 for j in range(i): 

267 d = get_distance(events[i], events[j], metric) 

268 simmat[i, j] = d 

269 simmat[j, i] = d 

270 

271 return simmat 

272 

273 

274def load_similarity_matrix(fname): 

275 ''' 

276 Load a binary similarity matrix from file 

277 ''' 

278 

279 simmat = num.load(fname) 

280 return simmat