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-26 16:25 +0000

1import math 

2import numpy as num 

3from pyrocko import moment_tensor, orthodrome 

4from grond.meta import GrondError 

5 

6 

7def get_distance_mt_l2(eventi, eventj): 

8 ''' 

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

10 ''' 

11 

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 

16 

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 

21 

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 

25 

26 d = 0.5 * math.sqrt(d) 

27 

28 return d 

29 

30 

31def get_distance_mt_l1(eventi, eventj): 

32 ''' 

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

34 ''' 

35 

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 

45 

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) 

49 

50 d = 0.5 * math.sqrt(d) 

51 

52 return d 

53 

54 

55def get_distance_mt_cos(eventi, eventj): 

56 ''' 

57 Inner product among two moment tensors. 

58 

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

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

61 ''' 

62 

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) 

70 

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) 

78 

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 

87 

88 if innerproduct >= 1.0: 

89 innerproduct = 1.0 

90 elif innerproduct <= -1.0: 

91 innerproduct = -1.0 

92 

93 d = 0.5 * (1 - innerproduct) 

94 

95 return d 

96 

97 

98def get_distance_mt_weighted_cos(eventi, eventj, ws): 

99 ''' 

100 Weighted moment tensor distance. 

101 

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) 

111 

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) 

119 

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 

128 

129 if innerproduct >= 1.0: 

130 innerproduct = 1.0 

131 

132 elif innerproduct <= -1.0: 

133 innerproduct = -1.0 

134 

135 d = 0.5 * (1.0 - innerproduct) 

136 

137 return d 

138 

139 

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 ''' 

144 

145 mti = eventi.moment_tensor 

146 mtj = eventj.moment_tensor 

147 

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

149 if d > 1.: 

150 d = 1. 

151 

152 return d 

153 

154 

155def get_distance_hypo(eventi, eventj): 

156 ''' 

157 Normalized Euclidean hypocentral distance, assuming flat earth to combine 

158 epicentral distance and depth difference. 

159 

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 

165 

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

167 

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) 

173 

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) 

178 

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

180 

181 d = hypo_distance_km / maxdist_km 

182 if d >= 1.: 

183 d = 1. 

184 

185 return d 

186 

187 

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. 

193 

194 a_lats, a_lons, b_lats, b_lons = \ 

195 eventi.north, eventi.east, eventj.north, eventj.east 

196 

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

198 

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) 

204 

205 distance_km = distance_m / 1000. 

206 

207 d = distance_km / maxdist_km 

208 if d >= 1.: 

209 d = 1. 

210 

211 return d 

212 

213 

214def get_distance_mt_triangle_diagram(eventi, eventj): 

215 ''' 

216 Scalar product among principal axes (?). 

217 ''' 

218 

219 mti = eventi.moment_tensor 

220 mtj = eventj.moment_tensor 

221 

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 

231 

232 if dotprod >= 1.: 

233 dotprod == 1. 

234 

235 d = 1. - dotprod 

236 return d 

237 

238 

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} 

249 

250 

251metrics = sorted(metric_funcs.keys()) 

252 

253 

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 ''' 

259 

260 try: 

261 func = metric_funcs[metric] 

262 except KeyError: 

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

264 

265 return func(eventi, eventj) 

266 

267 

268def compute_similarity_matrix(events, metric): 

269 ''' 

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

271 the desired metric 

272 

273 :param events: list of pyrocko events 

274 :param metric: metric type (string) 

275 

276 :returns: similarity matrix as NumPy array 

277 ''' 

278 

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 

286 

287 return simmat 

288 

289 

290def load_similarity_matrix(fname): 

291 ''' 

292 Load a binary similarity matrix from file 

293 ''' 

294 

295 simmat = num.load(fname) 

296 return simmat