Coverage for /usr/local/lib/python3.11/dist-packages/grond/scenario.py: 90%

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

2# 

3# The Grond Developers, 21st Century 

4import logging 

5import os 

6import shutil 

7import os.path as op 

8 

9from pyrocko import gf, scenario, util, guts 

10from pyrocko.gui import marker as pmarker 

11import math 

12 

13import grond 

14 

15 

16DEFAULT_STATIC_STORE = 'ak135_static' 

17DEFAULT_WAVEFORM_STORE = 'global2s' 

18 

19logger = logging.getLogger('grond.scenario') 

20km = 1e3 

21 

22 

23class GrondScenario(object): 

24 

25 def __init__(self, project_dir, 

26 center_lat=23., center_lon=52., radius=230*km, 

27 problem=None, observations=None): 

28 

29 if observations is None: 

30 observations = [] 

31 

32 self.project_dir = project_dir 

33 self.data_dir = op.join('data', 'scenario') 

34 

35 self.center_lat = center_lat 

36 self.center_lon = center_lon 

37 self.radius = radius 

38 

39 self.problem = problem 

40 self.observations = observations 

41 

42 self.rebuild = False 

43 self._scenario = None 

44 

45 def get_gf_stores_dir(self): 

46 return op.join(self.project_dir, 'gf_stores') 

47 

48 def create_project_dir(self, force=False): 

49 prj_dir = self.project_dir 

50 

51 if op.exists(prj_dir) and not force: 

52 raise grond.GrondError( 

53 'Directory "%s" already exists! Use --force to overwrite' 

54 % prj_dir) 

55 elif op.exists(prj_dir) and force: 

56 logger.info('Overwriting directory %s.', prj_dir) 

57 shutil.rmtree(prj_dir) 

58 

59 util.ensuredir(prj_dir) 

60 

61 @property 

62 def stores_wanted(self): 

63 target_generators = self.get_scenario().target_generators 

64 return set([t.store_id for t in target_generators]) 

65 

66 def symlink_gfstores(self, engine): 

67 logger.info('Symlinking Green\'s function stores...') 

68 

69 for store_id in self.stores_wanted: 

70 store = engine.get_store(store_id) 

71 dtarget = op.join(self.get_gf_stores_dir(), store_id) 

72 if not op.exists(dtarget): 

73 os.symlink(store.store_dir, dtarget) 

74 

75 def add_observation(self, observation): 

76 logger.info('Adding %s.' % observation.__class__.__name__) 

77 self.observations.append(observation) 

78 

79 def set_problem(self, problem): 

80 self.problem = problem 

81 

82 def get_dataset_config(self): 

83 events_path = op.join(self.data_dir, 'events.txt') 

84 

85 dataset_config = grond.DatasetConfig( 

86 events_path=events_path) 

87 

88 for obs in self.observations: 

89 obs.update_dataset_config(dataset_config, self.data_dir) 

90 return dataset_config 

91 

92 def set_config(self, filename): 

93 self.config_fn = filename 

94 

95 def get_scenario(self): 

96 if not self._scenario: 

97 if self.rebuild: 

98 scenario_file = op.join(self.project_dir, 'data', 'scenario', 

99 'scenario.yml') 

100 self._scenario = guts.load(filename=scenario_file) 

101 

102 else: 

103 if not self.observations: 

104 raise AttributeError('No observations set,' 

105 ' use .add_observation(Observation)') 

106 if not self.problem: 

107 raise AttributeError('No Source Problem set,' 

108 ' use .set_problem(Problem).') 

109 

110 self._scenario = scenario.ScenarioGenerator( 

111 center_lat=self.center_lat, 

112 center_lon=self.center_lon, 

113 radius=self.radius, 

114 target_generators=[obs.get_scenario_target_generator() 

115 for obs in self.observations], 

116 source_generator= self.problem.get_scenario_source_generator()) # noqa 

117 

118 return self._scenario 

119 

120 def create_scenario( 

121 self, 

122 force=False, 

123 interactive=True, 

124 gf_store_superdirs=None, 

125 make_map=True): 

126 

127 logger.info('Creating scenario...') 

128 

129 scenario = self.get_scenario() 

130 self.create_project_dir(force) 

131 util.ensuredir(self.get_gf_stores_dir()) 

132 

133 data_dir = op.join(self.project_dir, self.data_dir) 

134 util.ensuredir(data_dir) 

135 

136 scenario.dump(filename=op.join(data_dir, 'scenario.yml')) 

137 

138 if gf_store_superdirs is None: 

139 engine1 = gf.LocalEngine( 

140 use_config=True, 

141 store_superdirs=[self.get_gf_stores_dir()]) 

142 else: 

143 engine1 = gf.LocalEngine( 

144 use_config=False, 

145 store_superdirs=gf_store_superdirs) 

146 

147 scenario.init_modelling(engine=engine1) 

148 

149 scenario.ensure_gfstores( 

150 interactive=interactive) 

151 self.symlink_gfstores(engine1) 

152 

153 engine2 = gf.LocalEngine( 

154 use_config=False, 

155 store_superdirs=[self.get_gf_stores_dir()]) 

156 

157 scenario.init_modelling(engine=engine2) 

158 

159 scenario.dump_data(path=data_dir) 

160 if make_map: 

161 scenario.make_map(op.join(self.project_dir, 'scenario_map.pdf')) 

162 

163 shutil.move(op.join(data_dir, 'sources.yml'), 

164 op.join(data_dir, 'scenario_sources.yml')) 

165 

166 markers = scenario.get_onsets() 

167 marker_path = op.join(data_dir, 'picks', 'picks.markers') 

168 if markers: 

169 util.ensuredirs(marker_path) 

170 pmarker.save_markers(markers, marker_path) 

171 

172 def get_grond_config(self): 

173 engine_config = grond.EngineConfig( 

174 gf_stores_from_pyrocko_config=False, 

175 gf_store_superdirs=['gf_stores']) 

176 

177 optimiser_config = grond.HighScoreOptimiserConfig() 

178 

179 config = grond.Config( 

180 rundir_template=op.join('runs', '${problem_name}.grun'), 

181 dataset_config=self.get_dataset_config(), 

182 target_groups=[obs.get_grond_target_group() 

183 for obs in self.observations], 

184 problem_config=self.problem.get_grond_problem_config(), 

185 optimiser_config=optimiser_config, 

186 engine_config=engine_config) 

187 

188 config.set_basepath(self.project_dir) 

189 return config 

190 

191 def get_grond_config_path(self): 

192 return op.join('config', 'scenario.gronf') 

193 

194 def create_grond_files(self): 

195 logger.info('Creating Grond configuration for %s.' 

196 % ' and '.join([obs.name for obs in self.observations])) 

197 

198 config_path = op.join(self.project_dir, self.get_grond_config_path()) 

199 util.ensuredirs(config_path) 

200 grond.write_config(self.get_grond_config(), config_path) 

201 

202 def build( 

203 self, 

204 force=False, 

205 interactive=False, 

206 gf_store_superdirs=None, 

207 make_map=True): 

208 

209 logger.info('Building scenario...') 

210 

211 self.create_scenario( 

212 force=force, 

213 interactive=interactive, 

214 gf_store_superdirs=gf_store_superdirs, 

215 make_map=make_map) 

216 

217 self.create_grond_files() 

218 

219 

220class Observation(object): 

221 

222 name = 'Observation Prototype' 

223 

224 def __init__(self, store_id, *args, **kwargs): 

225 self.store_id = store_id 

226 

227 def update_dataset_config(self, dataset_config, data_dir): 

228 return dataset_config 

229 

230 def get_scenario_target_generator(self): 

231 pass 

232 

233 def get_grond_target_group(self): 

234 pass 

235 

236 

237class WaveformObservation(Observation): 

238 

239 name = 'seismic waveforms' 

240 

241 def __init__(self, store_id=DEFAULT_WAVEFORM_STORE, nstations=25, 

242 stationxml_paths=None, stations_paths=None): 

243 self.nstations = nstations 

244 self.store_id = store_id 

245 self.stationxml_paths = stationxml_paths 

246 self.stations_paths = stations_paths 

247 

248 def update_dataset_config(self, dataset_config, data_dir): 

249 ds = dataset_config 

250 ds.waveform_paths = [op.join(data_dir, 'waveforms')] 

251 ds.stations_stationxml_paths = [ 

252 op.join(data_dir, 'meta', 'waveform_response.xml')] 

253 ds.responses_stationxml_paths = [ 

254 op.join(data_dir, 'meta', 'waveform_response.xml')] 

255 return ds 

256 

257 def get_scenario_target_generator(self): 

258 if self.stationxml_paths or self.stations_paths: 

259 station_generator = scenario.station.ImportStationGenerator( 

260 stations_paths=self.stations_paths, 

261 stations_stationxml_paths=self.stationxml_paths) 

262 else: 

263 station_generator = scenario.station.RandomStationGenerator( 

264 nstations=self.nstations) 

265 return scenario.targets.WaveformGenerator( 

266 station_generators=[station_generator], 

267 store_id=self.store_id, 

268 tabulated_phases_from_store=True, 

269 tabulated_phases_noise_scale=0.3, 

270 seismogram_quantity='displacement') 

271 

272 def get_grond_target_group(self): 

273 return grond.WaveformTargetGroup( 

274 normalisation_family='time_domain', 

275 path='all', 

276 distance_min=10*km, 

277 distance_max=1000*km, 

278 channels=['Z', 'R', 'T'], 

279 interpolation='multilinear', 

280 store_id=self.store_id, 

281 misfit_config=grond.WaveformMisfitConfig( 

282 tmin='vel_surface:5.5', 

283 tmax='vel_surface:3.0', 

284 fmin=0.01, 

285 fmax=0.1)) 

286 

287 

288class InSARObservation(Observation): 

289 

290 name = 'InSAR displacement' 

291 

292 def __init__(self, store_id=DEFAULT_STATIC_STORE): 

293 self.store_id = store_id 

294 

295 def update_dataset_config(self, dataset_config, data_dir): 

296 dataset_config.kite_scene_paths = [op.join(data_dir, 'insar')] 

297 return dataset_config 

298 

299 def get_scenario_target_generator(self): 

300 return scenario.targets.InSARGenerator( 

301 mask_water=False, 

302 store_id=self.store_id) 

303 

304 def get_grond_target_group(self): 

305 logger.info('''Hint: 

306 

307Use kite's `spool` to configure the InSAR scene's quadtree and covariance! 

308''') 

309 return grond.SatelliteTargetGroup( 

310 normalisation_family='insar_target', 

311 path='all', 

312 interpolation='multilinear', 

313 store_id=self.store_id, 

314 kite_scenes=['*all'], 

315 misfit_config=grond.SatelliteMisfitConfig( 

316 optimise_orbital_ramp=True, 

317 ranges={ 

318 'offset': '-0.5 .. 0.5', 

319 'ramp_north': '-1e-4 .. 1e-4', 

320 'ramp_east': '-1e-4 .. 1e-4' 

321 }) 

322 ) 

323 

324 

325class GNSSCampaignObservation(Observation): 

326 

327 name = 'GNSS campaign' 

328 

329 def __init__(self, store_id=DEFAULT_STATIC_STORE, nstations=25): 

330 self.store_id = store_id 

331 self.nstations = nstations 

332 

333 def update_dataset_config(self, dataset_config, data_dir): 

334 dataset_config.gnss_campaign_paths = [op.join(data_dir, 'gnss')] 

335 return dataset_config 

336 

337 def get_scenario_target_generator(self): 

338 return scenario.targets.GNSSCampaignGenerator( 

339 station_generators=[scenario.targets.RandomStationGenerator( 

340 nstations=self.nstations)], 

341 store_id=self.store_id) 

342 

343 def get_grond_target_group(self): 

344 return grond.GNSSCampaignTargetGroup( 

345 gnss_campaigns=['*all'], 

346 normalisation_family='gnss_target', 

347 path='all', 

348 interpolation='multilinear', 

349 store_id=self.store_id, 

350 misfit_config=grond.GNSSCampaignMisfitConfig() 

351 ) 

352 

353 

354class SourceProblem(object): 

355 

356 def __init__(self, magnitude_min=6., magnitude_max=7., radius=10*km, 

357 nevents=1): 

358 self.magnitude_min = magnitude_min 

359 self.magnitude_max = magnitude_max 

360 self.radius = radius 

361 self.nevents = nevents 

362 

363 def get_scenario_source_generator(self): 

364 return 

365 

366 def get_grond_problem_config(self): 

367 return 

368 

369 

370class DCSourceProblem(SourceProblem): 

371 

372 def get_scenario_source_generator(self): 

373 return scenario.sources.DCSourceGenerator( 

374 magnitude_min=self.magnitude_min, 

375 magnitude_max=self.magnitude_max, 

376 radius=self.radius, 

377 depth_min=5*km, 

378 depth_max=20*km, 

379 nevents=self.nevents) 

380 

381 def get_grond_problem_config(self): 

382 s2 = math.sqrt(2) 

383 return grond.CMTProblemConfig( 

384 name_template='cmt_${event_name}', 

385 distance_min=2.*km, 

386 mt_type='full', 

387 ranges=dict( 

388 time=gf.Range(-5.0, 5.0, relative='add'), 

389 north_shift=gf.Range(-15*km, 15*km), 

390 east_shift=gf.Range(-15*km, 15*km), 

391 depth=gf.Range(5*km, 20*km), 

392 magnitude=gf.Range(self.magnitude_min, self.magnitude_max), 

393 rmnn=gf.Range(-1.0, 1.0), 

394 rmee=gf.Range(-1.0, 1.0), 

395 rmdd=gf.Range(-1.0, 1.0), 

396 rmne=gf.Range(-s2, s2), 

397 rmnd=gf.Range(-s2, s2), 

398 rmed=gf.Range(-s2, s2), 

399 duration=gf.Range(1.0, 15.0)) 

400 ) 

401 

402 

403class RectangularSourceProblem(SourceProblem): 

404 

405 def get_scenario_source_generator(self): 

406 return scenario.sources.RectangularSourceGenerator( 

407 magnitude_min=self.magnitude_min, 

408 magnitude_max=self.magnitude_max, 

409 depth_min=5*km, 

410 depth_max=15*km, 

411 radius=self.radius, 

412 nevents=self.nevents) 

413 

414 def get_grond_problem_config(self): 

415 return grond.RectangularProblemConfig( 

416 name_template='rect_source_${event_name}', 

417 decimation_factor=1, 

418 ranges=dict( 

419 north_shift=gf.Range(-15*km, 15*km), 

420 east_shift=gf.Range(-15*km, 15*km), 

421 depth=gf.Range(5*km, 15*km), 

422 length=gf.Range(20*km, 40*km), 

423 width=gf.Range(5*km, 10*km), 

424 dip=gf.Range(0, 90), 

425 strike=gf.Range(0, 180), 

426 rake=gf.Range(0, 90), 

427 slip=gf.Range(.5, 10), 

428 nucleation_x=gf.Range(-1., 1.), 

429 nucleation_y=gf.Range(-1., 1.), 

430 time=gf.Range(-5.0, 5.0, relative='add')) 

431 )