Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/scenario/targets/waveform.py: 92%

239 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2024-03-07 11:54 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Synthetic seismic waveform generator. 

8''' 

9 

10import hashlib 

11import math 

12import logging 

13import numpy as num 

14 

15from os import path as op 

16from functools import reduce 

17 

18from pyrocko.guts import Float, List, Bool, Dict, String 

19from pyrocko.gui.snuffler.marker import PhaseMarker, EventMarker 

20from pyrocko import gf, model, util, trace, io 

21from pyrocko.io.io_common import FileSaveError 

22from pyrocko import pile 

23 

24from ..station import StationGenerator, RandomStationGenerator 

25from .base import TargetGenerator, NoiseGenerator 

26from ..error import ScenarioError 

27 

28 

29DEFAULT_STORE_ID = 'global_2s' 

30 

31logger = logging.getLogger('pyrocko.scenario.targets.waveform') 

32guts_prefix = 'pf.scenario' 

33 

34 

35class WaveformNoiseGenerator(NoiseGenerator): 

36 

37 def get_time_increment(self, deltat): 

38 return deltat * 1024 

39 

40 def get_intersecting_snippets(self, deltat, codes, tmin, tmax): 

41 raise NotImplementedError() 

42 

43 def add_noise(self, tr): 

44 for ntr in self.get_intersecting_snippets( 

45 tr.deltat, tr.nslc_id, tr.tmin, tr.tmax): 

46 tr.add(ntr) 

47 

48 

49class WhiteNoiseGenerator(WaveformNoiseGenerator): 

50 

51 scale = Float.T(default=1e-6) 

52 

53 def get_seed_offset2(self, deltat, iw, codes): 

54 m = hashlib.sha1(('%e %i %s.%s.%s.%s' % ((deltat, iw) + codes)) 

55 .encode('utf8')) 

56 return int(m.hexdigest(), base=16) % 10000000 

57 

58 def get_intersecting_snippets(self, deltat, codes, tmin, tmax): 

59 tinc = self.get_time_increment(deltat) 

60 iwmin = int(math.floor(tmin / tinc)) 

61 iwmax = int(math.floor(tmax / tinc)) 

62 

63 trs = [] 

64 for iw in range(iwmin, iwmax+1): 

65 seed_offset = self.get_seed_offset2(deltat, iw, codes) 

66 rstate = self.get_rstate(seed_offset) 

67 

68 n = int(round(tinc // deltat)) 

69 

70 trs.append(trace.Trace( 

71 codes[0], codes[1], codes[2], codes[3], 

72 deltat=deltat, 

73 tmin=iw*tinc, 

74 ydata=rstate.normal(loc=0, scale=self.scale, size=n))) 

75 

76 return trs 

77 

78 

79class WaveformGenerator(TargetGenerator): 

80 

81 station_generators = List.T( 

82 StationGenerator.T(), 

83 default=[RandomStationGenerator.D()], 

84 help='List of StationGenerators.') 

85 

86 noise_generator = WaveformNoiseGenerator.T( 

87 default=WhiteNoiseGenerator.D(), 

88 help='Add Synthetic noise on the waveforms.') 

89 

90 store_id = gf.StringID.T( 

91 default=DEFAULT_STORE_ID, 

92 help='The GF store to use for forward-calculations.') 

93 

94 seismogram_quantity = gf.QuantityType.T( 

95 default='displacement') 

96 

97 nsl_to_store_id = Dict.T( 

98 String.T(), gf.StringID.T(), 

99 help='Selectively use different GF stores for different stations.') 

100 

101 vmin_cut = Float.T( 

102 default=2000., 

103 help='Minimum velocity to seismic velocity to consider in the model.') 

104 vmax_cut = Float.T( 

105 default=8000., 

106 help='Maximum velocity to seismic velocity to consider in the model.') 

107 

108 fmin = Float.T( 

109 default=0.01, 

110 help='Minimum frequency/wavelength to resolve in the' 

111 ' synthetic waveforms.') 

112 

113 tabulated_phases = List.T( 

114 gf.meta.TPDef.T(), optional=True, 

115 help='Define seismic phases to be calculated.') 

116 

117 tabulated_phases_from_store = Bool.T( 

118 default=False, 

119 help='Calculate seismic phase arrivals for all travel-time tables ' 

120 'defined in GF store.') 

121 

122 tabulated_phases_noise_scale = Float.T( 

123 default=0.0, 

124 help='Standard deviation of normally distributed noise added to ' 

125 'calculated phase arrivals.') 

126 

127 taper = trace.Taper.T( 

128 optional=True, 

129 help='Time domain taper applied to synthetic waveforms.') 

130 

131 compensate_synthetic_offsets = Bool.T( 

132 default=False, 

133 help='Center synthetic trace amplitudes using mean of waveform tips.') 

134 

135 tinc = Float.T( 

136 optional=True, 

137 help='Time increment of waveforms.') 

138 

139 continuous = Bool.T( 

140 default=True, 

141 help='Only produce traces that intersect with events.') 

142 

143 def __init__(self, *args, **kwargs): 

144 super(WaveformGenerator, self).__init__(*args, **kwargs) 

145 self._targets = [] 

146 self._piles = {} 

147 

148 def _get_pile(self, path): 

149 apath = op.abspath(path) 

150 assert op.isdir(apath) 

151 

152 if apath not in self._piles: 

153 fns = util.select_files( 

154 [apath], show_progress=False) 

155 

156 p = pile.Pile() 

157 if fns: 

158 p.load_files(fns, fileformat='mseed', show_progress=False) 

159 

160 self._piles[apath] = p 

161 

162 return self._piles[apath] 

163 

164 def get_stations(self): 

165 stations = [] 

166 for station_generator in self.station_generators: 

167 stations.extend(station_generator.get_stations()) 

168 return stations 

169 

170 def get_distance_range(self, sources): 

171 distances = num.array( 

172 [sg.get_distance_range(sources) 

173 for sg in self.station_generators]) 

174 return (distances[:, 0].min(), distances[:, 1].max()) 

175 

176 def get_targets(self): 

177 if self._targets: 

178 return self._targets 

179 

180 for station in self.get_stations(): 

181 channel_data = [] 

182 channels = station.get_channels() 

183 if channels: 

184 for channel in channels: 

185 channel_data.append([ 

186 channel.name, 

187 channel.azimuth, 

188 channel.dip]) 

189 

190 else: 

191 for c_name in ['BHZ', 'BHE', 'BHN']: 

192 channel_data.append([ 

193 c_name, 

194 model.guess_azimuth_from_name(c_name), 

195 model.guess_dip_from_name(c_name)]) 

196 

197 nsl = (station.network, station.station, station.location) 

198 for c_name, c_azi, c_dip in channel_data: 

199 target = gf.Target( 

200 codes=nsl + (c_name,), 

201 quantity=self.seismogram_quantity, 

202 lat=station.lat, 

203 lon=station.lon, 

204 north_shift=station.north_shift, 

205 east_shift=station.east_shift, 

206 depth=station.depth, 

207 store_id=self.nsl_to_store_id.get( 

208 '.'.join(nsl), self.store_id), 

209 optimization='enable', 

210 interpolation='nearest_neighbor', 

211 azimuth=c_azi, 

212 dip=c_dip) 

213 

214 self._targets.append(target) 

215 

216 return self._targets 

217 

218 def get_time_range(self, sources): 

219 dmin, dmax = self.get_distance_range(sources) 

220 

221 times = num.array([source.time for source in sources], 

222 dtype=util.get_time_dtype()) 

223 

224 tmin_events = num.min(times) 

225 tmax_events = num.max(times) 

226 

227 tmin = tmin_events + dmin / self.vmax_cut - 10.0 / self.fmin 

228 tmax = tmax_events + dmax / self.vmin_cut + 10.0 / self.fmin 

229 

230 return tmin, tmax 

231 

232 def get_codes_to_deltat(self, engine, sources): 

233 deltats = {} 

234 

235 targets = self.get_targets() 

236 for source in sources: 

237 for target in targets: 

238 deltats[target.codes] = engine.get_store( 

239 target.store_id).config.deltat 

240 

241 return deltats 

242 

243 def get_useful_time_increment(self, engine, sources): 

244 _, dmax = self.get_distance_range(sources) 

245 tinc = dmax / self.vmin_cut + 2.0 / self.fmin 

246 

247 deltats = set(self.get_codes_to_deltat(engine, sources).values()) 

248 deltat = reduce(util.lcm, deltats) 

249 tinc = int(round(tinc / deltat)) * deltat 

250 return tinc 

251 

252 def get_relevant_sources(self, sources, tmin, tmax): 

253 dmin, dmax = self.get_distance_range(sources) 

254 trange = tmax - tmin 

255 tmax_pad = trange + tmax + dmin / self.vmax_cut 

256 tmin_pad = tmin - (dmax / self.vmin_cut + trange) 

257 

258 return [s for s in sources if s.time < tmax_pad and s.time > tmin_pad] 

259 

260 def get_waveforms(self, engine, sources, tmin, tmax): 

261 

262 sources_relevant = self.get_relevant_sources(sources, tmin, tmax) 

263 if not (self.continuous or sources_relevant): 

264 return [] 

265 

266 trs = {} 

267 tts = util.time_to_str 

268 

269 for nslc, deltat in self.get_codes_to_deltat(engine, sources).items(): 

270 tr_tmin = round(tmin / deltat) * deltat 

271 tr_tmax = (round(tmax / deltat)-1) * deltat 

272 nsamples = int(round((tr_tmax - tr_tmin) / deltat)) + 1 

273 

274 tr = trace.Trace( 

275 *nslc, 

276 tmin=tr_tmin, 

277 ydata=num.zeros(nsamples), 

278 deltat=deltat) 

279 

280 self.noise_generator.add_noise(tr) 

281 

282 trs[nslc] = tr 

283 

284 logger.debug('Forward modelling waveforms between %s - %s...' 

285 % (tts(tmin, format='%Y-%m-%d_%H-%M-%S'), 

286 tts(tmax, format='%Y-%m-%d_%H-%M-%S'))) 

287 

288 if not sources_relevant: 

289 return list(trs.values()) 

290 

291 targets = self.get_targets() 

292 response = engine.process(sources_relevant, targets) 

293 for source, target, res in response.iter_results( 

294 get='results'): 

295 

296 if isinstance(res, gf.SeismosizerError): 

297 logger.warning( 

298 'Out of bounds! \nTarget: %s\nSource: %s\n' % ( 

299 '.'.join(target.codes)), source) 

300 continue 

301 

302 tr = res.trace.pyrocko_trace() 

303 

304 candidate = trs[target.codes] 

305 if not candidate.overlaps(tr.tmin, tr.tmax): 

306 continue 

307 

308 if self.compensate_synthetic_offsets: 

309 tr.ydata -= (num.mean(tr.ydata[-3:-1]) + 

310 num.mean(tr.ydata[1:3])) / 2. 

311 

312 if self.taper: 

313 tr.taper(self.taper) 

314 

315 candidate.add(tr) 

316 trs[target.codes] = candidate 

317 

318 return list(trs.values()) 

319 

320 def get_onsets(self, engine, sources, *args, **kwargs): 

321 

322 targets = {t.codes[:3]: t for t in self.get_targets()} 

323 

324 markers = [] 

325 for source in sources: 

326 ev = source.pyrocko_event() 

327 markers.append(EventMarker(ev)) 

328 for nsl, target in targets.items(): 

329 store = engine.get_store(target.store_id) 

330 if self.tabulated_phases: 

331 tabulated_phases = self.tabulated_phases 

332 

333 elif self.tabulated_phases_from_store: 

334 tabulated_phases = store.config.tabulated_phases 

335 else: 

336 tabulated_phases = [] 

337 

338 for phase in tabulated_phases: 

339 t = store.t(phase.id, source, target) 

340 if not t: 

341 continue 

342 

343 noise_scale = self.tabulated_phases_noise_scale 

344 if noise_scale != 0.0: 

345 t += num.random.normal(scale=noise_scale) 

346 

347 t += source.time 

348 markers.append( 

349 PhaseMarker( 

350 phasename=phase.id, 

351 tmin=t, 

352 tmax=t, 

353 event=source.pyrocko_event(), 

354 nslc_ids=(nsl+('*',),) 

355 ) 

356 ) 

357 return markers 

358 

359 def ensure_data(self, engine, sources, path, tmin=None, tmax=None): 

360 self.ensure_waveforms(engine, sources, path, tmin, tmax) 

361 self.ensure_responses(path) 

362 

363 def ensure_waveforms(self, engine, sources, path, tmin=None, tmax=None): 

364 

365 path_waveforms = op.join(path, 'waveforms') 

366 util.ensuredir(path_waveforms) 

367 

368 p = self._get_pile(path_waveforms) 

369 

370 nslc_ids = set(target.codes for target in self.get_targets()) 

371 

372 def have_waveforms(tmin, tmax): 

373 trs_have = p.all( 

374 tmin=tmin, tmax=tmax, 

375 load_data=False, degap=False, 

376 trace_selector=lambda tr: tr.nslc_id in nslc_ids) 

377 

378 return any(tr.data_len() > 0 for tr in trs_have) 

379 

380 def add_files(paths): 

381 p.load_files(paths, fileformat='mseed', show_progress=False) 

382 

383 path_traces = op.join( 

384 path_waveforms, 

385 '%(wmin_year)s', 

386 '%(wmin_month)s', 

387 '%(wmin_day)s', 

388 'waveform_%(network)s_%(station)s_' + 

389 '%(location)s_%(channel)s_%(tmin)s_%(tmax)s.mseed') 

390 

391 tmin_all, tmax_all = self.get_time_range(sources) 

392 tmin = tmin if tmin is not None else tmin_all 

393 tmax = tmax if tmax is not None else tmax_all 

394 tts = util.time_to_str 

395 

396 tinc = self.tinc or self.get_useful_time_increment(engine, sources) 

397 tmin = num.floor(tmin / tinc) * tinc 

398 tmax = num.ceil(tmax / tinc) * tinc 

399 

400 nwin = int(round((tmax - tmin) / tinc)) 

401 

402 pbar = None 

403 try: 

404 for iwin in range(nwin): 

405 tmin_win = tmin + iwin*tinc 

406 tmax_win = tmin + (iwin+1)*tinc 

407 

408 if have_waveforms(tmin_win, tmax_win): 

409 continue 

410 

411 if pbar is None: 

412 pbar = util.progressbar( 

413 'Generating waveforms', (nwin-iwin)) 

414 

415 pbar.update(iwin) 

416 

417 trs = self.get_waveforms(engine, sources, tmin_win, tmax_win) 

418 

419 try: 

420 wpaths = io.save( 

421 trs, path_traces, 

422 additional=dict( 

423 wmin_year=tts(tmin_win, format='%Y'), 

424 wmin_month=tts(tmin_win, format='%m'), 

425 wmin_day=tts(tmin_win, format='%d'), 

426 wmin=tts(tmin_win, format='%Y-%m-%d_%H-%M-%S'), 

427 wmax_year=tts(tmax_win, format='%Y'), 

428 wmax_month=tts(tmax_win, format='%m'), 

429 wmax_day=tts(tmax_win, format='%d'), 

430 wmax=tts(tmax_win, format='%Y-%m-%d_%H-%M-%S'))) 

431 

432 for wpath in wpaths: 

433 logger.debug('Generated file: %s' % wpath) 

434 

435 add_files(wpaths) 

436 

437 except FileSaveError as e: 

438 raise ScenarioError(str(e)) 

439 

440 finally: 

441 if pbar is not None: 

442 pbar.finish() 

443 

444 def ensure_responses(self, path): 

445 from pyrocko.io import stationxml 

446 

447 path_responses = op.join(path, 'meta') 

448 util.ensuredir(path_responses) 

449 

450 fn_stationxml = op.join(path_responses, 'waveform_response.xml') 

451 i = 0 

452 while op.exists(fn_stationxml): 

453 fn_stationxml = op.join( 

454 path_responses, 'waveform_response-%s.xml' % i) 

455 i += 1 

456 

457 logger.debug('Writing waveform meta information to StationXML...') 

458 

459 stations = self.get_stations() 

460 sxml = stationxml.FDSNStationXML.from_pyrocko_stations(stations) 

461 

462 sunit = stationxml.quantity_to_units[self.seismogram_quantity] 

463 

464 response = stationxml.Response( 

465 instrument_sensitivity=stationxml.Sensitivity( 

466 value=1., 

467 frequency=1., 

468 input_units=stationxml.Units(sunit), 

469 output_units=stationxml.Units('COUNT')), 

470 stage_list=[]) 

471 

472 for net, station, channel in sxml.iter_network_station_channels(): 

473 channel.response = response 

474 

475 sxml.dump_xml(filename=fn_stationxml) 

476 

477 def add_map_artists(self, engine, sources, automap): 

478 automap.add_stations(self.get_stations())