Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/waveform_phase_ratio/target.py: 44%

1import logging 

2 

3import numpy as num 

4 

5from pyrocko.guts import Float, Tuple, String, Bool 

6from pyrocko import gf 

7 

8from ..base import ( 

9 MisfitTarget, TargetGroup, MisfitResult) 

10from . import measure as fm 

11from grond import dataset 

12from grond.meta import has_get_plot_classes 

13 

14from ..waveform.target import StoreIDSelector 

15 

16guts_prefix = 'grond' 

17logger = logging.getLogger('grond.targets.waveform_phase_ratio.target') 

18 

19 

20def log_exclude(target, reason): 

21 logger.debug('Excluding potential target %s: %s' % ( 

22 target.string_id(), reason)) 

23 

24 

25class PhaseRatioTargetGroup(TargetGroup): 

26 

27 ''' 

28 Generate targets to compare ratios or log ratios of two seismic phases. 

29 

30 misfit = | a_obs / (a_obs + b_obs) - a_syn / (a_syn + b_syn) | 

31 

32 or 

33 

34 misfit = | log(a_obs / (a_obs + b_obs) + waterlevel) - 

35 log(a_syn / (a_syn + b_syn) + waterlevel) | 

36 ''' 

37 

38 distance_min = Float.T(optional=True) 

39 distance_max = Float.T(optional=True) 

40 distance_3d_min = Float.T(optional=True) 

41 distance_3d_max = Float.T(optional=True) 

42 depth_min = Float.T(optional=True) 

43 depth_max = Float.T(optional=True) 

44 measure_a = fm.FeatureMeasure.T() 

45 measure_b = fm.FeatureMeasure.T() 

46 interpolation = gf.InterpolationMethod.T() 

47 store_id = gf.StringID.T(optional=True) 

48 store_id_selector = StoreIDSelector.T( 

49 optional=True, 

50 help='select GF store based on event-station geometry.') 

51 

52 fit_log_ratio = Bool.T( 

53 default=True, 

54 help='If true, compare synthetic and observed log ratios') 

55 

56 fit_log_ratio_waterlevel = Float.T( 

57 default=0.01, 

58 help='Waterlevel added to both ratios when comparing on logarithmic ' 

59 'scale, to avoid log(0)') 

60 

61 def get_targets(self, ds, event, default_path='none'): 

62 logger.debug('Selecting phase ratio targets...') 

63 origin = event 

64 targets = [] 

65 

66 for st in ds.get_stations(): 

67 blacklisted = False 

68 for measure in [self.measure_a, self.measure_b]: 

69 for cha in measure.channels: 

70 if ds.is_blacklisted((st.nsl() + (cha,))): 

71 blacklisted = True 

72 

73 if self.store_id_selector: 

74 store_id = self.store_id_selector.get_store_id( 

75 event, st, cha) 

76 else: 

77 store_id = self.store_id 

78 

79 target = PhaseRatioTarget( 

80 codes=st.nsl(), 

81 lat=st.lat, 

82 lon=st.lon, 

83 north_shift=st.north_shift, 

84 east_shift=st.east_shift, 

85 depth=st.depth, 

86 interpolation=self.interpolation, 

87 store_id=store_id, 

88 measure_a=self.measure_a, 

89 measure_b=self.measure_b, 

90 manual_weight=self.weight, 

91 normalisation_family=self.normalisation_family, 

92 path=self.path or default_path, 

93 backazimuth=0.0, 

94 fit_log_ratio=self.fit_log_ratio, 

95 fit_log_ratio_waterlevel=self.fit_log_ratio_waterlevel) 

96 

97 if blacklisted: 

98 log_exclude(target, 'blacklisted') 

99 continue 

100 

101 if self.distance_min is not None and \ 

102 target.distance_to(origin) < self.distance_min: 

103 log_exclude(target, 'distance < distance_min') 

104 continue 

105 

106 if self.distance_max is not None and \ 

107 target.distance_to(origin) > self.distance_max: 

108 log_exclude(target, 'distance > distance_max') 

109 continue 

110 

111 if self.distance_3d_min is not None and \ 

112 target.distance_3d_to(origin) < self.distance_3d_min: 

113 log_exclude(target, 'distance_3d < distance_3d_min') 

114 continue 

115 

116 if self.distance_3d_max is not None and \ 

117 target.distance_3d_to(origin) > self.distance_3d_max: 

118 log_exclude(target, 'distance_3d > distance_3d_max') 

119 continue 

120 

121 if self.depth_min is not None and \ 

122 target.depth < self.depth_min: 

123 log_exclude(target, 'depth < depth_min') 

124 continue 

125 

126 if self.depth_max is not None and \ 

127 target.depth > self.depth_max: 

128 log_exclude(target, 'depth > depth_max') 

129 continue 

130 

131 bazi, _ = target.azibazi_to(origin) 

132 target.backazimuth = bazi 

133 target.set_dataset(ds) 

134 targets.append(target) 

135 

136 return targets 

137 

138 

139class PhaseRatioResult(MisfitResult): 

140 a_obs = Float.T(optional=True) 

141 b_obs = Float.T(optional=True) 

142 a_syn = Float.T(optional=True) 

143 b_syn = Float.T(optional=True) 

144 

145 

146@has_get_plot_classes 

147class PhaseRatioTarget(gf.Location, MisfitTarget): 

148 

149 ''' 

150 Target to compare ratios or log ratios of two seismic phases. 

151 

152 misfit = | a_obs / (a_obs + b_obs) - a_syn / (a_syn + b_syn) | 

153 

154 ''' 

155 

156 codes = Tuple.T( 

157 3, String.T(), 

158 help='network, station, location codes.') 

159 

160 store_id = gf.StringID.T( 

161 help='ID of Green\'s function store to use for the computation.') 

162 

163 backazimuth = Float.T(optional=True) 

164 

165 interpolation = gf.InterpolationMethod.T() 

166 

167 measure_a = fm.FeatureMeasure.T() 

168 measure_b = fm.FeatureMeasure.T() 

169 

170 fit_log_ratio = Bool.T( 

171 default=True, 

172 help='if true, compare synthetic and observed log ratios') 

173 

174 fit_log_ratio_waterlevel = Float.T( 

175 default=0.01, 

176 help='Waterlevel added to both ratios when comparing on logarithmic ' 

177 'scale, to avoid log(0)') 

178 

179 can_bootstrap_weights = True 

180 

181 def __init__(self, **kwargs): 

182 gf.Location.__init__(self, **kwargs) 

183 MisfitTarget.__init__(self, **kwargs) 

184 

185 @classmethod 

186 def get_plot_classes(cls): 

187 from . import plot 

188 plots = super(PhaseRatioTarget, cls).get_plot_classes() 

189 plots.extend(plot.get_plot_classes()) 

190 return plots 

191 

192 def string_id(self): 

193 return '.'.join(x for x in (self.path,) + self.codes) 

194 

195 def get_plain_targets(self, engine, source): 

196 return self.prepare_modelling(engine, source, None) 

197 

198 def prepare_modelling(self, engine, source, targets): 

199 modelling_targets = [] 

200 for measure in [self.measure_a, self.measure_b]: 

201 modelling_targets.extend(measure.get_modelling_targets( 

202 self.codes, self.lat, self.lon, self.depth, self.store_id, 

203 self.backazimuth)) 

204 

205 return modelling_targets 

206 

207 def finalize_modelling( 

208 self, engine, source, modelling_targets, modelling_results): 

209 

210 ds = self.get_dataset() 

211 

212 try: 

213 imt = 0 

214 amps = [] 

215 for measure in [self.measure_a, self.measure_b]: 

216 nmt_this = measure.get_nmodelling_targets() 

217 amp_obs, _ = measure.evaluate( 

218 engine, source, 

219 modelling_targets[imt:imt+nmt_this], 

220 dataset=ds) 

221 

222 amp_syn, _ = measure.evaluate( 

223 engine, source, 

224 modelling_targets[imt:imt+nmt_this], 

225 trs=[r.trace.pyrocko_trace() 

226 for r 

227 in modelling_results[imt:imt+nmt_this]]) 

228 

229 amps.append((amp_obs, amp_syn)) 

230 

231 imt += nmt_this 

232 

233 (a_obs, a_syn), (b_obs, b_syn) = amps 

234 

235 eps = self.fit_log_ratio_waterlevel 

236 if self.fit_log_ratio: 

237 res_a = num.log(a_obs / (a_obs + b_obs) + eps) \ 

238 - num.log(a_syn / (a_syn + b_syn) + eps) 

239 else: 

240 res_a = a_obs / (a_obs + b_obs) - a_syn / (a_syn + b_syn) 

241 

242 misfit = num.abs(res_a) 

243 norm = 1.0 

244 

245 result = PhaseRatioResult( 

246 misfits=num.array([[misfit, norm]], dtype=float), 

247 a_obs=a_obs, 

248 b_obs=b_obs, 

249 a_syn=a_syn, 

250 b_syn=b_syn) 

251 

252 return result 

253 

254 except dataset.NotFound as e: 

255 logger.debug(str(e)) 

256 return gf.SeismosizerError('No waveform data: %s' % str(e)) 

257 

258 

259__all__ = ''' 

260 PhaseRatioTargetGroup 

261 PhaseRatioTarget 

262 PhaseRatioResult 

263'''.split()