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

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

2# 

3# The Grond Developers, 21st Century 

4import logging 

5 

6import numpy as num 

7 

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

9from pyrocko import gf 

10 

11from ..base import ( 

12 MisfitTarget, TargetGroup, MisfitResult) 

13from . import measure as fm 

14from grond import dataset 

15from grond.meta import has_get_plot_classes 

16 

17from ..waveform.target import StoreIDSelector 

18 

19guts_prefix = 'grond' 

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

21 

22 

23def log_exclude(target, reason): 

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

25 target.string_id(), reason)) 

26 

27 

28class PhaseRatioTargetGroup(TargetGroup): 

29 

30 ''' 

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

32 

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

34 

35 or 

36 

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

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

39 ''' 

40 

41 distance_min = Float.T(optional=True) 

42 distance_max = Float.T(optional=True) 

43 distance_3d_min = Float.T(optional=True) 

44 distance_3d_max = Float.T(optional=True) 

45 depth_min = Float.T(optional=True) 

46 depth_max = Float.T(optional=True) 

47 measure_a = fm.FeatureMeasure.T() 

48 measure_b = fm.FeatureMeasure.T() 

49 interpolation = gf.InterpolationMethod.T() 

50 store_id_selector = StoreIDSelector.T( 

51 optional=True, 

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

53 

54 fit_log_ratio = Bool.T( 

55 default=True, 

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

57 

58 fit_log_ratio_waterlevel = Float.T( 

59 default=0.01, 

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

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

62 

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

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

65 origin = event 

66 targets = [] 

67 

68 for st in ds.get_stations(): 

69 excluded = False 

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

71 for cha in measure.channels: 

72 if ds.is_excluded((st.nsl() + (cha,))): 

73 excluded = True 

74 

75 if self.store_id_selector: 

76 store_id = self.store_id_selector.get_store_id( 

77 event, st, cha) 

78 else: 

79 store_id = self.store_id 

80 

81 target = PhaseRatioTarget( 

82 codes=st.nsl(), 

83 lat=st.lat, 

84 lon=st.lon, 

85 north_shift=st.north_shift, 

86 east_shift=st.east_shift, 

87 depth=st.depth, 

88 interpolation=self.interpolation, 

89 store_id=store_id, 

90 measure_a=self.measure_a, 

91 measure_b=self.measure_b, 

92 manual_weight=self.weight, 

93 normalisation_family=self.normalisation_family, 

94 path=self.path or default_path, 

95 backazimuth=0.0, 

96 fit_log_ratio=self.fit_log_ratio, 

97 fit_log_ratio_waterlevel=self.fit_log_ratio_waterlevel) 

98 

99 if excluded: 

100 log_exclude(target, 'excluded') 

101 continue 

102 

103 if self.distance_min is not None and \ 

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

105 log_exclude(target, 'distance < distance_min') 

106 continue 

107 

108 if self.distance_max is not None and \ 

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

110 log_exclude(target, 'distance > distance_max') 

111 continue 

112 

113 if self.distance_3d_min is not None and \ 

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

115 log_exclude(target, 'distance_3d < distance_3d_min') 

116 continue 

117 

118 if self.distance_3d_max is not None and \ 

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

120 log_exclude(target, 'distance_3d > distance_3d_max') 

121 continue 

122 

123 if self.depth_min is not None and \ 

124 target.depth < self.depth_min: 

125 log_exclude(target, 'depth < depth_min') 

126 continue 

127 

128 if self.depth_max is not None and \ 

129 target.depth > self.depth_max: 

130 log_exclude(target, 'depth > depth_max') 

131 continue 

132 

133 bazi, _ = target.azibazi_to(origin) 

134 target.backazimuth = bazi 

135 target.set_dataset(ds) 

136 targets.append(target) 

137 

138 return targets 

139 

140 def is_gf_store_appropriate(self, store, depth_range): 

141 ok = TargetGroup.is_gf_store_appropriate( 

142 self, store, depth_range) 

143 ok &= self._is_gf_store_appropriate_check_extent( 

144 store, depth_range) 

145 ok &= self._is_gf_store_appropriate_check_sample_rate( 

146 store, depth_range) 

147 return ok 

148 

149 

150class PhaseRatioResult(MisfitResult): 

151 a_obs = Float.T(optional=True) 

152 b_obs = Float.T(optional=True) 

153 a_syn = Float.T(optional=True) 

154 b_syn = Float.T(optional=True) 

155 

156 

157@has_get_plot_classes 

158class PhaseRatioTarget(gf.Location, MisfitTarget): 

159 

160 ''' 

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

162 

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

164 

165 ''' 

166 

167 codes = Tuple.T( 

168 3, String.T(), 

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

170 

171 store_id = gf.StringID.T( 

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

173 

174 backazimuth = Float.T(optional=True) 

175 

176 interpolation = gf.InterpolationMethod.T() 

177 

178 measure_a = fm.FeatureMeasure.T() 

179 measure_b = fm.FeatureMeasure.T() 

180 

181 fit_log_ratio = Bool.T( 

182 default=True, 

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

184 

185 fit_log_ratio_waterlevel = Float.T( 

186 default=0.01, 

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

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

189 

190 can_bootstrap_weights = True 

191 

192 def __init__(self, **kwargs): 

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

194 MisfitTarget.__init__(self, **kwargs) 

195 

196 @classmethod 

197 def get_plot_classes(cls): 

198 from . import plot 

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

200 plots.extend(plot.get_plot_classes()) 

201 return plots 

202 

203 def string_id(self): 

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

205 

206 def get_plain_targets(self, engine, source): 

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

208 

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

210 modelling_targets = [] 

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

212 modelling_targets.extend(measure.get_modelling_targets( 

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

214 self.backazimuth)) 

215 

216 return modelling_targets 

217 

218 def finalize_modelling( 

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

220 

221 ds = self.get_dataset() 

222 

223 try: 

224 imt = 0 

225 amps = [] 

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

227 nmt_this = measure.get_nmodelling_targets() 

228 amp_obs, _ = measure.evaluate( 

229 engine, source, 

230 modelling_targets[imt:imt+nmt_this], 

231 dataset=ds) 

232 

233 amp_syn, _ = measure.evaluate( 

234 engine, source, 

235 modelling_targets[imt:imt+nmt_this], 

236 trs=[r.trace.pyrocko_trace() 

237 for r 

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

239 

240 amps.append((amp_obs, amp_syn)) 

241 

242 imt += nmt_this 

243 

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

245 

246 eps = self.fit_log_ratio_waterlevel 

247 if self.fit_log_ratio: 

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

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

250 else: 

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

252 

253 misfit = num.abs(res_a) 

254 norm = 1.0 

255 

256 result = PhaseRatioResult( 

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

258 a_obs=a_obs, 

259 b_obs=b_obs, 

260 a_syn=a_syn, 

261 b_syn=b_syn) 

262 

263 return result 

264 

265 except dataset.NotFound as e: 

266 logger.debug(str(e)) 

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

268 

269 

270__all__ = ''' 

271 PhaseRatioTargetGroup 

272 PhaseRatioTarget 

273 PhaseRatioResult 

274'''.split()