Coverage for /usr/local/lib/python3.11/dist-packages/grond/scenario.py: 90%
187 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-11-27 15:15 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2024-11-27 15:15 +0000
1# https://pyrocko.org/grond - GPLv3
2#
3# The Grond Developers, 21st Century
4import logging
5import os
6import shutil
7import os.path as op
9from pyrocko import gf, scenario, util, guts
10from pyrocko.gui import marker as pmarker
11import math
13import grond
16DEFAULT_STATIC_STORE = 'ak135_static'
17DEFAULT_WAVEFORM_STORE = 'global2s'
19logger = logging.getLogger('grond.scenario')
20km = 1e3
23class GrondScenario(object):
25 def __init__(self, project_dir,
26 center_lat=23., center_lon=52., radius=230*km,
27 problem=None, observations=None):
29 if observations is None:
30 observations = []
32 self.project_dir = project_dir
33 self.data_dir = op.join('data', 'scenario')
35 self.center_lat = center_lat
36 self.center_lon = center_lon
37 self.radius = radius
39 self.problem = problem
40 self.observations = observations
42 self.rebuild = False
43 self._scenario = None
45 def get_gf_stores_dir(self):
46 return op.join(self.project_dir, 'gf_stores')
48 def create_project_dir(self, force=False):
49 prj_dir = self.project_dir
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)
59 util.ensuredir(prj_dir)
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])
66 def symlink_gfstores(self, engine):
67 logger.info('Symlinking Green\'s function stores...')
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)
75 def add_observation(self, observation):
76 logger.info('Adding %s.' % observation.__class__.__name__)
77 self.observations.append(observation)
79 def set_problem(self, problem):
80 self.problem = problem
82 def get_dataset_config(self):
83 events_path = op.join(self.data_dir, 'events.txt')
85 dataset_config = grond.DatasetConfig(
86 events_path=events_path)
88 for obs in self.observations:
89 obs.update_dataset_config(dataset_config, self.data_dir)
90 return dataset_config
92 def set_config(self, filename):
93 self.config_fn = filename
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)
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).')
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
118 return self._scenario
120 def create_scenario(
121 self,
122 force=False,
123 interactive=True,
124 gf_store_superdirs=None,
125 make_map=True):
127 logger.info('Creating scenario...')
129 scenario = self.get_scenario()
130 self.create_project_dir(force)
131 util.ensuredir(self.get_gf_stores_dir())
133 data_dir = op.join(self.project_dir, self.data_dir)
134 util.ensuredir(data_dir)
136 scenario.dump(filename=op.join(data_dir, 'scenario.yml'))
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)
147 scenario.init_modelling(engine=engine1)
149 scenario.ensure_gfstores(
150 interactive=interactive)
151 self.symlink_gfstores(engine1)
153 engine2 = gf.LocalEngine(
154 use_config=False,
155 store_superdirs=[self.get_gf_stores_dir()])
157 scenario.init_modelling(engine=engine2)
159 scenario.dump_data(path=data_dir)
160 if make_map:
161 scenario.make_map(op.join(self.project_dir, 'scenario_map.pdf'))
163 shutil.move(op.join(data_dir, 'sources.yml'),
164 op.join(data_dir, 'scenario_sources.yml'))
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)
172 def get_grond_config(self):
173 engine_config = grond.EngineConfig(
174 gf_stores_from_pyrocko_config=False,
175 gf_store_superdirs=['gf_stores'])
177 optimiser_config = grond.HighScoreOptimiserConfig()
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)
188 config.set_basepath(self.project_dir)
189 return config
191 def get_grond_config_path(self):
192 return op.join('config', 'scenario.gronf')
194 def create_grond_files(self):
195 logger.info('Creating Grond configuration for %s.'
196 % ' and '.join([obs.name for obs in self.observations]))
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)
202 def build(
203 self,
204 force=False,
205 interactive=False,
206 gf_store_superdirs=None,
207 make_map=True):
209 logger.info('Building scenario...')
211 self.create_scenario(
212 force=force,
213 interactive=interactive,
214 gf_store_superdirs=gf_store_superdirs,
215 make_map=make_map)
217 self.create_grond_files()
220class Observation(object):
222 name = 'Observation Prototype'
224 def __init__(self, store_id, *args, **kwargs):
225 self.store_id = store_id
227 def update_dataset_config(self, dataset_config, data_dir):
228 return dataset_config
230 def get_scenario_target_generator(self):
231 pass
233 def get_grond_target_group(self):
234 pass
237class WaveformObservation(Observation):
239 name = 'seismic waveforms'
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
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
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')
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))
288class InSARObservation(Observation):
290 name = 'InSAR displacement'
292 def __init__(self, store_id=DEFAULT_STATIC_STORE):
293 self.store_id = store_id
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
299 def get_scenario_target_generator(self):
300 return scenario.targets.InSARGenerator(
301 mask_water=False,
302 store_id=self.store_id)
304 def get_grond_target_group(self):
305 logger.info('''Hint:
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 )
325class GNSSCampaignObservation(Observation):
327 name = 'GNSS campaign'
329 def __init__(self, store_id=DEFAULT_STATIC_STORE, nstations=25):
330 self.store_id = store_id
331 self.nstations = nstations
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
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)
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 )
354class SourceProblem(object):
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
363 def get_scenario_source_generator(self):
364 return
366 def get_grond_problem_config(self):
367 return
370class DCSourceProblem(SourceProblem):
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)
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 )
403class RectangularSourceProblem(SourceProblem):
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)
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 )