Browse Source

Merge pull request 'Update of source inversion to use the mole architecture' (#6) from si_update into master

Reviewed-on: #6
master
Malte Metz 2 weeks ago
parent
commit
8355a84f8e
  1. 7
      src/apps/ids.py
  2. 19
      src/apps/mole.py
  3. 139
      src/data/examples/grond_dynamic_source.gronf
  4. 173
      src/data/examples/grond_mt_source.gronf
  5. 4
      src/gm/config.py
  6. 4
      src/ls/config.py
  7. 84
      src/plot/ids.py
  8. 292
      src/si/config.py
  9. 6
      src/si/ids/core.py
  10. 185
      src/si/invert.py
  11. 103
      src/si/module.py
  12. 36
      src/store.py
  13. 9
      test/si/test_ids.py

7
src/apps/ids.py

@ -120,7 +120,12 @@ def command_init(args):
if args.rundir is None:
raise IDSError('running directoy needs to be given')
IDSRunner.init(rundir=args.rundir, force=args.force)
fns = IDSRunner.init(rundir=args.rundir, force=args.force)
logger.info('(1) configure settings in config file e.g.:\n %s'
% '\n '.join(fns))
logger.info('(2) run "ids go <config_file>" in directory: "%s"'
% args.rundir)
def command_go(args):

19
src/apps/mole.py

@ -152,16 +152,17 @@ def command_go(args):
parser.add_argument(
'--offline', dest='offline', action='store_true',
help='If given, offline not realtime mode is done')
help='If given, offline not realtime mode is done. '
'default is %(default)s')
parser.add_argument(
'--playback_sim', dest='playback', action='store_true',
help='If given, offline data is replayed as in realtime')
help='If given, offline data is replayed as in realtime. '
'default is %(default)s')
parser.add_argument(
'--event', dest='event_id', action='store',
help='If given, evaluate local data in ./data/<event_id>. '
'default is %(default)s')
help='If given, evaluate local data in ./data/<event_id>.')
parser.add_argument(
'--data', dest='data_dir', action='store',
@ -189,9 +190,13 @@ def command_go(args):
'--ls_conf', dest='ls_conf', action='store',
help='LandSlide config file, if config/ls.conf shall not be used')
parser.add_argument(
'--create-report', dest='create_report', action='store_true',
help='If given, Grond and IDS reports are created.')
parser, args = cl_parse('go', args, setup=setup)
storage = store.ResultStorage(store_dir='.')
storage = store.ResultStorage(store_dir=op.abspath('.'))
storage.open(
conf_dir=args.config,
si_conf=args.si_conf,
@ -209,7 +214,9 @@ def command_go(args):
playback=args.playback,
configs=storage.configs)
problem.run(force=args.force)
problem.set_parent(storage)
problem.run(force=args.force, create_report=args.create_report)
# TODO capture logs, errors etc

139
src/data/examples/grond_dynamic_source.gronf

@ -0,0 +1,139 @@
%YAML 1.1
--- !grond.Config
# Is automatically updated by mole
path_prefix: /local/home
rundir_template: runs_fastmodel_syn/${problem_name}.grun
# Dataset config is automatically update by mole
dataset_config: !grond.DatasetConfig
stations_stationxml_paths: ['path/to/stationxml.xml']
events_path: 'path/to/event/file.yaml'
waveform_paths: ['path/to/waveforms/directory']
responses_stationxml_paths: ['path/to/responses.xml']
extend_incomplete: true
picks_paths: []
blacklist_paths: ['path/to/blacklist.txt']
blacklist: [TU.AYVA, KO.BODT, KO.DATO, KO.TVSB, KO.CAVI]
whitelist_paths: []
target_groups:
- !grond.WaveformTargetGroup
normalisation_family: td
path: td.p
weight: 1.0
interpolation: nearest_neighbor
store_id: ak135_2000km_1Hz
distance_min: 0.0
distance_max: 150000.0
include: [NO.IKAR, NO.094A, NO.SAMU]
exclude: []
channels: [Z, R]
misfit_config: !grond.WaveformMisfitConfig
quantity: displacement
fmin: 0.01
fmax: 0.05
ffactor: 1.5
tmin: '0.'
tmax: '80.'
tfade: 20.0
domain: time_domain
norm_exponent: 1
tautoshift_max: 5.0
autoshift_penalty_max: 0.3
- !grond.WaveformTargetGroup
normalisation_family: td
path: td.pacc
weight: 1.0
interpolation: nearest_neighbor
store_id: ak135_2000km_1Hz
distance_min: 0.0
distance_max: 150000.0
include: [KO.GMLD, KO.KUSD, TK.0911, TK.0918, TK.0920, TK.3511, TK.3512, TK.3514,
TK.3517, TK.3524, TK.3523, TK.3528, TK.3533, TK.3538]
exclude: []
channels: [Z, R]
misfit_config: !grond.WaveformMisfitConfig
quantity: acceleration
fmin: 0.01
fmax: 0.05
ffactor: 1.5
tmin: '0.'
tmax: '80.'
tfade: 20.0
domain: time_domain
norm_exponent: 1
tautoshift_max: 5.0
autoshift_penalty_max: 0.3
- !grond.WaveformTargetGroup
normalisation_family: td
path: td.pbb
weight: 1.0
interpolation: nearest_neighbor
store_id: ak135_500km_1Hz
distance_min: 100000.0
distance_max: 200000.0
include: [HL.APE, HL.PRK, HL.PRK, KO.MLSB, KO.YAYO, KO.DAT, KO.BODT, KO.DKL,
KO.YER, KO.GORD, KO.NAZL, KO.SOMA, KO.YAYO, TU.CAM, TU.TURN,
TU.BDRM, TU.SUL, TU.STEP, TU.KOCA, TU.MULA, TU.BUHA, TU.AKHS,
TU.ZEDA, TU.AYVA, TU.MANT, TU.KTT, TU.IZMD, TU.ESEN, TU.IZMR,
TU.YAZI, TU.KIRA]
exclude: []
channels: [Z, R]
misfit_config: !grond.WaveformMisfitConfig
quantity: displacement
fmin: 0.01
fmax: 0.05
ffactor: 1.5
tmin: '0.'
tmax: '80.'
tfade: 20.0
domain: time_domain
norm_exponent: 1
tautoshift_max: 5.0
autoshift_penalty_max: 0.3
problem_config: !grond.DynamicRuptureProblemConfig
name_template: 'dynamic_rupture_run'
norm_exponent: 1
ranges:
depth: 1000 .. 10000
dip: 0 .. 90
east_shift: -25000 .. 25000
gamma: 0.7 .. 0.9
length: 10000 .. 60000
north_shift: -25000 .. 25000
nucleation_x: -1 .. 1
nucleation_y: -1 .. 1
nx: 9 .. 9
ny: 5 .. 5
rake: -180 .. 180
slip: 0.1 .. 7
strike: -180 .. 180
time: -15 .. 15 | add
width: 5000 .. 20000
anchor: top
decimation_factor: 1
distance_min: 1000.0
nthreads: 4
adaptive_resolution: 'off'
adaptive_start: 0
pure_shear: true
smooth_rupture: true
point_source_target_balancing: true
analyser_configs:
- !grond.TargetBalancingAnalyserConfig
niterations: 1000
use_reference_magnitude: false
optimiser_config: !grond.HighScoreOptimiserConfig
sampler_phases: []
# - !grond.InjectionSamplerPhase
# niterations: 2236
# ntries_preconstrain_limit: 1000
# sources_path: /data/local/home/pyrocko/projects/samos_2020/grond/Samos_20201030_115127/sources_inject_test0.yaml
chain_length_factor: 8.0
nbootstrap: 20
engine_config: !grond.EngineConfig
gf_stores_from_pyrocko_config: false
gf_store_superdirs: [/absolute/path/to/pyrocko/gf/stores]
gf_store_dirs: []
event_names: [event_name]
event_names_exclude: []

173
src/data/examples/grond_mt_source.gronf

@ -0,0 +1,173 @@
%YAML 1.1
--- !grond.Config
# Is automatically updated by mole
path_prefix: /local/home
rundir_template: runs_fastmodel_syn/${problem_name}.grun
# Is automatically updated by mole
dataset_config: !grond.DatasetConfig
stations_stationxml_paths: ['path/to/stationxml.xml']
events_path: 'path/to/event/file.yaml'
waveform_paths: ['path/to/waveforms/directory']
responses_stationxml_paths: ['path/to/responses.xml']
extend_incomplete: true
picks_paths: []
blacklist_paths: ['path/to/blacklist.txt']
blacklist: [TU.AYVA, KO.BODT, KO.DATO, KO.TVSB, KO.CAVI]
whitelist_paths: []
target_groups:
- !grond.WaveformTargetGroup
normalisation_family: td
path: td.rayleighacc
weight: 0.75
interpolation: nearest_neighbor
store_id: ak135_2000km_1Hz
distance_min: 0.0
distance_max: 100000.0
include: [KO.GMLD, KO.KUSD, KO.CESE, KO.DIDI, TK.0911, TK.0918, TK.0920, TK.3511,
TK.3512, TK.3514, TK.3517, TK.3524, TK.3523, TK.3528, TK.3533, TK.3536, TK.3538]
exclude: []
channels: [R, Z]
misfit_config: !grond.WaveformMisfitConfig
quantity: acceleration
fmin: 0.01
fmax: 0.03
ffactor: 1.5
tmin: '{stored:any_P}'
tmax: '{vel_surface:2.0}'
domain: time_domain
norm_exponent: 1
tautoshift_max: 2.5
autoshift_penalty_max: 0.05
- !grond.WaveformTargetGroup
normalisation_family: td
path: td.loveacc
weight: 0.75
interpolation: nearest_neighbor
store_id: ak135_2000km_1Hz
distance_min: 0.0
distance_max: 100000.0
include: [KO.GMLD, KO.KUSD, KO.CESE, KO.DIDI, TK.0911, TK.0918, TK.0920, TK.3511,
TK.3512, TK.3514, TK.3517, TK.3524, TK.3523, TK.3528, TK.3533, TK.3536, TK.3538]
exclude: []
channels: [T]
misfit_config: !grond.WaveformMisfitConfig
quantity: acceleration
fmin: 0.01
fmax: 0.03
ffactor: 1.5
tmin: '{stored:any_P}'
tmax: '{vel_surface:2.0}'
domain: time_domain
norm_exponent: 1
tautoshift_max: 2.5
autoshift_penalty_max: 0.05
- !grond.WaveformTargetGroup
normalisation_family: td
path: td.rayleigh
weight: 1.0
interpolation: nearest_neighbor
store_id: ak135_2000km_1Hz
distance_min: 100000.0
distance_max: 900000.0
exclude: [TK.3539, TK.1016, TK.1027, TK.3525, TK.3512, TK.4814, TK.1019, TK.4823,
TK.0905, TK.0910, TK.3511, TK.4812, TK.4502, TK.1719, TK.1716, TK.4801, TK.4507,
TK.0913, TK.3526, TK.1017, TK.4509, TK.4511, TK.3521, TK.1028, TK.3519, TK.3503,
TK.4503, TK.3523, TK.3518, TK.0917, TK.0921, TK.3537, TK.0911, TK.0918, TK.1013,
TK.3506, TK.1704, TK.0916, TK.0919, TK.1724, TK.3514, TK.4815, TK.4819, TK.0920,
TK.3513, TK.1005, TK.4818, TK.0922, TK.4508, TK.3508, TK.4821, TK.1022, TK.3534,
TK.3517, TK.3524, TK.3516, TK.4506, TK.4822, TK.4809, TK.1008, TK.4810, TK.3522,
TK.3538, TK.3528, TK.1707, TK.1021, TK.3536, TK.3520, TK.4807, TK.1720, TK.3527,
TK.1015, TK.0912, TK.4808, TK.1026, TK.4501, TK.3533, TK.4817, TK.4806, KO.AKS,
KO.CANM, KO.CESE, KO.DATC, KO.DIDI, KO.EZNE, KO.FOCM, KO.GMLD, KO.GOMA, KO.KRBN,
KO.KUSD, KO.YENI, KO.YKAV, TU.BOZC, TU.CNKL, TU.IZMD, TU.IZMR]
channels: [R, Z]
misfit_config: !grond.WaveformMisfitConfig
quantity: displacement
fmin: 0.01
fmax: 0.03
ffactor: 1.5
tmin: '{stored:any_P}'
tmax: '{vel_surface:2.0}'
domain: time_domain
norm_exponent: 1
tautoshift_max: 2.5
autoshift_penalty_max: 0.05
- !grond.WaveformTargetGroup
normalisation_family: td
path: td.love
weight: 1.0
interpolation: nearest_neighbor
store_id: ak135_2000km_1Hz
distance_min: 100000.0
distance_max: 900000.0
exclude: [TK.3539, TK.1016, TK.1027, TK.3525, TK.3512, TK.4814, TK.1019, TK.4823,
TK.0905, TK.0910, TK.3511, TK.4812, TK.4502, TK.1719, TK.1716, TK.4801, TK.4507,
TK.0913, TK.3526, TK.1017, TK.4509, TK.4511, TK.3521, TK.1028, TK.3519, TK.3503,
TK.4503, TK.3523, TK.3518, TK.0917, TK.0921, TK.3537, TK.0911, TK.0918, TK.1013,
TK.3506, TK.1704, TK.0916, TK.0919, TK.1724, TK.3514, TK.4815, TK.4819, TK.0920,
TK.3513, TK.1005, TK.4818, TK.0922, TK.4508, TK.3508, TK.4821, TK.1022, TK.3534,
TK.3517, TK.3524, TK.3516, TK.4506, TK.4822, TK.4809, TK.1008, TK.4810, TK.3522,
TK.3538, TK.3528, TK.1707, TK.1021, TK.3536, TK.3520, TK.4807, TK.1720, TK.3527,
TK.1015, TK.0912, TK.4808, TK.1026, TK.4501, TK.3533, TK.4817, TK.4806, KO.AKS,
KO.CANM, KO.CESE, KO.DATC, KO.DIDI, KO.EZNE, KO.FOCM, KO.GMLD, KO.GOMA, KO.KRBN,
KO.KUSD, KO.YENI, KO.YKAV, TU.BOZC, TU.CNKL, TU.IZMD, TU.IZMR]
channels: [T]
misfit_config: !grond.WaveformMisfitConfig
quantity: displacement
fmin: 0.01
fmax: 0.033
ffactor: 1.5
tmin: '{stored:any_P}'
tmax: '{vel_surface:2.0}'
domain: time_domain
norm_exponent: 1
tautoshift_max: 2.5
autoshift_penalty_max: 0.05
problem_config: !grond.CMTProblemConfig
name_template: 'moment_tensor_run'
norm_exponent: 1
ranges:
depth: 2000 .. 30000
duration: 5 .. 30
east_shift: -25000 .. 25000
magnitude: 4 .. 8
north_shift: -25000 .. 25000
rmdd: -1.41421 .. 1.41421
rmed: -1 .. 1
rmee: -1.41421 .. 1.41421
rmnd: -1 .. 1
rmne: -1 .. 1
rmnn: -1.41421 .. 1.41421
time: -15 .. 15 | add
distance_min: 1000.0
mt_type: full
stf_type: HalfSinusoidSTF
nthreads: 1
analyser_configs:
- !grond.TargetBalancingAnalyserConfig
niterations: 1000
use_reference_magnitude: false
optimiser_config: !grond.HighScoreOptimiserConfig
sampler_phases: []
# - !grond.UniformSamplerPhase
# niterations: 1000
# ntries_preconstrain_limit: 1000
# - !grond.DirectedSamplerPhase
# niterations: 20000
# ntries_preconstrain_limit: 1000
# scatter_scale_begin: 2.0
# scatter_scale_end: 0.45
# starting_point: excentricity_compensated
# sampler_distribution: normal
# standard_deviation_estimator: median_density_single_chain
# ntries_sample_limit: 1000
chain_length_factor: 8.0
nbootstrap: 100
engine_config: !grond.EngineConfig
gf_stores_from_pyrocko_config: false
gf_store_superdirs: [/absolute/path/to/pyrocko/gf/stores]
gf_store_dirs: []
event_names: [Samos_20201030_115127]
event_names_exclude: []

4
src/gm/config.py

@ -13,5 +13,9 @@ def get_config(*args, **kwargs):
return GroundMotionConfig(*args, **kwargs)
def get_module_extra(self):
return []
class GroundMotionConfig(ModuleConfig):
pass

4
src/ls/config.py

@ -13,5 +13,9 @@ def get_config(*args, **kwargs):
return LandSlideConfig(*args, **kwargs)
def get_module_extra(self):
return []
class LandSlideConfig(ModuleConfig):
pass

84
src/plot/ids.py

@ -96,10 +96,15 @@ def _mplcmap_to_gmtcpt_code(mplcmap, steps=256):
Get gmt readable R/G/A code from a given matplotlib colormap.
:param mplcmap:
Name of the demanded matplotlib colormap
Name of the demanded matplotlib colormap.
:type mplcmap:
str
:param steps:
Number of discrete color steps.
:type steps:
optional, int
:returns:
Series of comma seperate R/G/B values for gmtpy usage.
:rtype:
@ -662,7 +667,8 @@ class RuptureMap(Map):
widths,
data,
filename,
interpolation='nearest'):
interpolation='nearest',
draw_zeros=True):
'''
Interpolate regular data onto topotile grid.
@ -694,6 +700,12 @@ class RuptureMap(Map):
``linear`` and ``nearest``.
:type interpolation:
optional, str
:param draw_zeros:
If ``True``, zeros are displayed as colored planes. If ``False``,
not. Default is ``True``.
:type draw_zeros:
optional, bool
'''
source = self.source
@ -717,6 +729,9 @@ class RuptureMap(Map):
t.data = interpolator(points_out).reshape(t.data.shape)
if not draw_zeros:
t.data[t.data == 0.] = num.nan
_save_grid(t.y(), t.x(), t.data, filename=filename)
def _irregular_data_to_grid(
@ -729,7 +744,8 @@ class RuptureMap(Map):
dips,
data,
filename,
interpolation='nearest'):
interpolation='nearest',
draw_zeros=True):
'''
Interpolate irregular data onto topotile grid.
@ -784,6 +800,12 @@ class RuptureMap(Map):
``linear`` and ``nearest``.
:type interpolation:
optional, str
:param draw_zeros:
If ``True``, zeros are displayed as colored planes. If ``False``,
not. Default is ``True``.
:type draw_zeros:
optional, bool
'''
from scipy.interpolate import griddata
@ -861,20 +883,23 @@ class RuptureMap(Map):
mask = num.ma.mask_or(mask, msk)
if interpolation == 'nearest':
data_interp[msk] = da[i]
# if interpolation == 'nearest':
# data_interp[msk] = da[i]
if interpolation == 'linear':
data_interp = griddata(
center_points,
da,
(meshgrid_lats, meshgrid_lons),
method=interpolation).flatten()
# if interpolation == 'linear':
data_interp = griddata(
center_points,
da,
(meshgrid_lats, meshgrid_lons),
method=interpolation).flatten()
data_interp[~mask] = num.nan
t.data = data_interp.reshape(t.data.shape)
if not draw_zeros:
t.data[t.data == 0.] = num.nan
_save_grid(t.y(), t.x(), t.data, filename=filename)
def patch_data_to_grid(self, data, *args, **kwargs):
@ -959,19 +984,22 @@ class RuptureMap(Map):
str
:param cbar:
If True, a colorbar corresponding to the grid data is added.
If ``True``, a colorbar corresponding to the grid data is added.
Keyword arguments are parsed to it.
:type cbar:
optional, bool
'''
kwargs_img = {}
kwargs_img['t'] = kwargs.pop('t', '0')
kwargs_img['C'] = kwargs.pop('C', 'my_%s.cpt' % cmap)
kwargs_img['E'] = kwargs.pop('E', '200')
kwargs_img['Q'] = kwargs.pop('Q', True)
kwargs_img['n'] = kwargs.pop('n', '+t0.0')
self.gmt.grdimage(
gridfile,
C='my_%s.cpt' % cmap,
E='200',
Q=True,
n='+t0.0',
*self.jxyr)
*self.jxyr,
**kwargs_img)
if cbar:
self.draw_colorbar(cmap=cmap, **kwargs)
@ -1060,6 +1088,7 @@ class RuptureMap(Map):
w_add = '+c'
kwargs['A'], kwargs['C'] = a_string, c_string
kwargs['S'] = kwargs.get('S', '10')
if style:
style = ',' + style
@ -1068,7 +1097,6 @@ class RuptureMap(Map):
self.gmt.grdcontour(
gridfile,
S='10',
W='a%gp,%s%s+s%s' % (pen_size*4, color, style, w_add),
*self.jxyr + args,
**kwargs)
@ -1119,7 +1147,6 @@ class RuptureMap(Map):
'D', 'j%s+w%gc/%gc+h+o%gc/%gc' % (a_str, w, h, dx, dy))
kwargs['F'] = kwargs.get(
'F', '+g238/236/230+c%g/%g/%g/%g' % (lgap, rgap, bgap, tgap))
print(kwargs['D'])
self.gmt.psscale(
*self.jxyr,
@ -1240,12 +1267,21 @@ class RuptureMap(Map):
cpt = [kwargs['cmap']]
tmp_grd_file = 'tmpdata.grd'
self.patch_data_to_grid(plot_data, tmp_grd_file)
kwargs_grd = {}
if 'draw_zeros' in kwargs:
kwargs_grd['draw_zeros'] = kwargs.get('draw_zeros', True)
del kwargs['draw_zeros']
if 'interpolation' in kwargs:
kwargs_grd['interpolation'] = kwargs.get('interpolation', True)
del kwargs['interpolation']
self.patch_data_to_grid(plot_data, tmp_grd_file, **kwargs_grd)
self.draw_image(tmp_grd_file, **kwargs)
clear_temp(gridfiles=[tmp_grd_file], cpts=cpt)
def draw_patch_parameter(self, attribute, **kwargs):
def draw_patch_parameter(self, attribute, draw_zeros=True, **kwargs):
'''
Draw an image of a chosen patch attribute.
@ -1268,7 +1304,7 @@ class RuptureMap(Map):
'label',
'%s [%s]' % (a, cbar_helper[a]['unit']))
self.draw_dynamic_data(data, **kwargs)
self.draw_dynamic_data(data, draw_zeros=draw_zeros, **kwargs)
def draw_time_contour(
self,
@ -1412,7 +1448,7 @@ class RuptureMap(Map):
self.draw_points(nlat, nlon, **kwargs)
def draw_dislocation(self, component='', **kwargs):
def draw_dislocation(self, component='', draw_zeros=True, **kwargs):
'''
Draw dislocation onto map at any time.
@ -1437,7 +1473,7 @@ class RuptureMap(Map):
kwargs['label'] = kwargs.get(
'label', 'u%s [m]' % (component))
self.draw_dynamic_data(data, **kwargs)
self.draw_dynamic_data(data, draw_zeros=draw_zeros, **kwargs)
def draw_dislocation_contour(
self,

292
src/si/config.py

@ -2,24 +2,17 @@
#
# The Developers, 21st Century
import logging
import numpy as num
import os.path as op
from pyrocko import model
from pyrocko.guts import Float, Tuple, String, Bool, List
from pyrocko.gf import Range
from pyrocko.gf.meta import QuantityType
from pyrocko.guts import Int
from grond import (
Path, read_config, Config,
DirectedSamplerPhase, HighScoreOptimiserConfig,
UniformSamplerPhase, InjectionSamplerPhase, EngineConfig,
DatasetConfig)
from grond.targets import WaveformMisfitConfig, WaveformTargetGroup
from grond.problems import CMTProblemConfig
from grond.analysers.target_balancing import TargetBalancingAnalyserConfig
from grond import Path, read_config, Config
from ewrica import util
from ewrica.module import ModuleConfig
from ewrica.dataset import tsumaps
from .invert import SourceInjectionConfig, GrondInverter, NEAMTHM18Limit
logger = logging.getLogger('ewrica.si.config')
@ -27,226 +20,105 @@ logger = logging.getLogger('ewrica.si.config')
km = 1e3
ps_config_fn = 'grond_mt_source.gronf'
ds_config_fn = 'grond_dynamic_source.gronf'
def get_config(*args, **kwargs):
return SourceInversionConfig(*args, **kwargs)
def get_module_extra():
return [util.data_file(
op.join('examples', fn)) for fn in (ps_config_fn, ds_config_fn)]
class InversionConfig(ModuleConfig):
config_path = Path.T(
grond_config_fn = Path.T(
default='',
help='Path to the grond inversion configuration file.')
config__ = Config.T(
grond_config__ = Config.T(
optional=True,
help='loaded grond inversion configuration object')
help='Loaded grond inversion configuration object.')
def __init__(self, *args, **kwargs):
ModuleConfig.__init__(self, *args, **kwargs)
self.config__ = None
injection_config = SourceInjectionConfig.T(
optional=True,
default=SourceInjectionConfig(),
help='Loaded injection configuration object.')
def generate_config(self, *args, **kwargs):
raise NotImplementedError
threads = Int.T(
default=1,
help='Number of threads.')
def get_config(self, *args, **kwargs):
if self.config__:
return self.config__
if self.config_path:
self.config__ = read_config(self.config_path)
return self.config__
self.generate_config(*args, **kwargs)
return self.config__
class PointSourceInversionConfig(InversionConfig):
distance_range = Tuple.T(
default=(100*km, 200*km))
fmin = Float.T(
default=0.01,
help='Minimum waveform frequency')
fmax = Float.T(
default=0.1,
help='Maximum waveform frequency')
quantity = QuantityType.T(
default='displacement',
help='Target quantity, where synthetic and real data are compared.')
mt_type = String.T(
default='full',
help=('Components of the full moment tensor, which are inverted:'
'"full", "deviatoric" or "dc".'))
inject_models_from_tsumaps = Bool.T(
default=True,
help='If true, models based on tsumaps probabilities are injected.')
store_id = String.T(
default='ak135_2000km_1Hz',
help='Greens function store id')
gf_store_superdirs = List.T(
Path.T(),
default=['.'],
help='Greens function store super directory')
event_file = Path.T(
default='',
help='Event file path')
event_id = String.T(
default='',
help='Event ID')
station_paths = List.T(
Path.T(),
default=[],
help='List of StationXML station and response files')
waveform_paths = List.T(
Path.T(),
default=[],
help='List of paths to waveform directories')
parent_path = Path.T(
default='.',
help='Path to parent directory. Other paths can be relative to it.')
blocklist_paths = List.T(
Path.T(),
parallel = Int.T(
default=1,
help='Number of parallel processes.')
base_path = Path.T(
optional=True,
default=None,
help='Path to txt file with blocked stations.')
def generate_config(self, *args, **kwargs):
target_groups = []
target_groups += [WaveformTargetGroup(
normalisation_family='td',
path='td.p',
distance_min=self.distance_range[0],
distance_max=self.distance_range[1],
channels=['Z', 'R'],
weight=1.0,
interpolation='nearest_neighbor',
store_id=self.store_id,
misfit_config=WaveformMisfitConfig(
quantity=self.quantity,
fmin=self.fmin,
fmax=self.fmax,
ffactor=1.5,
tmin='{stored:any_P}',
tmax='{stored:any_P}+50.',
domain='time_domain',
tautoshift_max=2.5,
autoshift_penalty_max=0.05,
norm_exponent=1))]
problem_conf = CMTProblemConfig(
name_template='',
ranges=dict(
time=Range(start=-15, stop=15, relative='add'),
north_shift=Range(start=-25.e3, stop=25.e3, step=1e3),
east_shift=Range(start=-25.e3, stop=25.e3, step=1e3),
depth=Range(start=1.e3, stop=25.e3, step=1e3),
magnitude=Range(start=4.0, stop=8.0, step=0.1),
rmnn=Range(start=-1.41421, stop=1.41421, step=0.05),
rmee=Range(start=-1.41421, stop=1.41421, step=0.05),
rmdd=Range(start=-1.41421, stop=1.41421, step=0.05),
rmne=Range(start=-1., stop=1., step=0.05),
rmnd=Range(start=-1., stop=1., step=0.05),
rmed=Range(start=-1., stop=1., step=0.05),
duration=Range(start=5., stop=20., step=1.)),
distance_min=1e3,
mt_type=self.mt_type,
norm_exponent=1)
anaylser_confs = []
anaylser_confs += [TargetBalancingAnalyserConfig(niterations=1000)]
optimiser_phases = []
optimiser_phases += [UniformSamplerPhase(niterations=1000)]
optimiser_phases += [DirectedSamplerPhase(
niterations=5000,
scatter_scale_begin=2.0,
scatter_scale_end=0.5)]
optimiser_conf = HighScoreOptimiserConfig(
nbootstrap=50, sampler_phases=optimiser_phases)
engine_conf = EngineConfig(
gf_stores_from_pyrocko_config=False,
gf_store_superdirs=self.gf_store_superdirs)
dataset_conf = DatasetConfig(
stations_stationxml_paths=self.station_paths,
responses_stationxml_paths=self.station_paths,
events_path=self.event_file,
waveform_paths=self.waveform_paths,
blacklist_paths=self.blocklist_paths or None,
apply_correction_factors=True,
extend_incomplete=True)
dataset_conf.set_basepath('/')
return Config(
rundir_template=op.join(
self.parent_path,
self.event_id, 'si', 'proc', '${problem_name}.grun'),
path_prefix='', # self.parent_path,
event_names=[self.event_id],
dataset_config=dataset_conf,
target_groups=target_groups,
problem_config=problem_conf,
analyser_configs=anaylser_confs,
optimiser_config=optimiser_conf,
engine_config=engine_conf)
def optimiser_conf_from_tsumaps(self, config):
optimiser_phases = []
optimiser_phases += [UniformSamplerPhase(niterations=1000)]
if self.inject_models_from_tsumaps:
events = model.load_events(self.event_file)
event = None
for ev in events:
if ev.name == self.event_id:
event = ev
if event is None:
raise ValueError(
'No event called %s found in the given event file' %
self.event_id)
xs_inject = tsumaps.tsumaps_to_models(
lat=event.lat, lon=event.lon,
ranges=config.problem_config.ranges,
n_models=1000, cum_prob_max=1.)
optimiser_phases += [InjectionSamplerPhase(
niterations=xs_inject.shape[0], xs_inject=xs_inject)]
optimiser_phases += [DirectedSamplerPhase(
niterations=5000,
scatter_scale_begin=2.0,
scatter_scale_end=0.5)]
config.optimiser_config = HighScoreOptimiserConfig(
nbootstrap=50, sampler_phases=optimiser_phases)
return config
help='Absolute path to main directory.')
def get_config(self, *args, **kwargs):
if self.config__ is not None:
return self.config__
report_path = Path.T(
optional=True,
help='If given, directory, where Grond reports are stored.')
def __init__(self, *args, **kwargs):
ModuleConfig.__init__(self, *args, **kwargs)
def load_grond_config_from_file(self):
if self.grond_config__ is None and self.grond_config_fn is not None:
self.grond_config = read_config(self.grond_config_fn)
if self.config_path:
config = read_config(self.config_path)
else:
config = self.generate_config(*args, **kwargs)
@property
def grond_config(self):
return self.grond_config__
if self.inject_models_from_tsumaps:
config = self.optimiser_conf_from_tsumaps(config)
@grond_config.setter
def grond_config(self, grond_config):
self.grond_config__ = grond_config
self.config__ = config
def get_config(self, *args, **kwargs):
return self
return config
def get_inverter(self, *args, **kwargs):
inverter = GrondInverter(
grond_config=self.grond_config,
injection_config=self.injection_config,
threads=self.threads,
parallel=self.parallel,
base_path=self.base_path,
report_path=self.report_path,
*args, **kwargs)
return inverter
class SourceInversionConfig(ModuleConfig):
pointsource = PointSourceInversionConfig.T(
default=PointSourceInversionConfig(),
pointsource_config = InversionConfig.T(
default=InversionConfig(grond_config_fn=ps_config_fn),
help='Configurations for point source inversion')
dynamic_rupture_config = InversionConfig.T(
default=InversionConfig(
grond_config_fn=ds_config_fn,
injection_config=SourceInjectionConfig(
asperity=['total', 'left', 'center', 'right'],
nucleation_x=num.array([-.66, 0., .66]),
nucleation_y=num.array([0., 0., 0.]),
tractions=num.array([0.5, 1., 2.5]) * 1e6,
neamthm18_limit=NEAMTHM18Limit(
quantity='nmodels_max',
target_value=7),
fault_database='EDSF',
length_factor=1.,
width_factor=1.,
include_dimension_uncertainty=False,
f_max=None,
v_min=1200.,
correct_event_duration=False)),
help='Configurations for dynamic rupture inversion')
# @property
# def paths(self):
# return dict(

6
src/si/ids/core.py

@ -187,11 +187,15 @@ class IDSRunner(Object):
os.makedirs(rundir)
os.makedirs(op.join(rundir, 'config'))
fns = [op.join(rundir, 'config', 'idsconfig_example.conf')]
runner = cls()
runner.config.dump(
filename=op.join(rundir, 'config', 'idsconfig_example.conf'),
filename=fns[0],
header=True)
return fns
def _prepare_ids_dir(self):
dirs = \
[self.waveformdir, self.gnssdir, self.insardir, self.faultdir] + \

185
src/si/invert.py

@ -32,7 +32,8 @@ d2r = num.pi / 180.
class NEAMTHM18Limit(Object):
'''Sets limiting factor and value for model extraction from NEAMTHM18 db
'''
Limiting factor and value for model extraction from NEAMTHM18 database.
'''
quantity = StringChoice.T(
choices=(
@ -42,6 +43,15 @@ class NEAMTHM18Limit(Object):
target_value = Float.T(
optional=True)
def get_target_value(self):
if self.target_value is None:
return None
if self.quantity == 'nmodels_max':
return int(self.target_value)
return self.target_value
class SourceInjectionConfig(Object):
nucleation_x = Array.T(
@ -82,7 +92,9 @@ class SourceInjectionConfig(Object):
optional=True)
neamthm18_limit = NEAMTHM18Limit.T(
default=NEAMTHM18Limit(quantity='min_gradient', target_value=0.01),
default=NEAMTHM18Limit(
quantity='probability_gradient_min',
target_value=0.01),
help='Sets the extracted models from NEAMTHM18 database.')
fault_database = StringChoice.T(
@ -112,7 +124,6 @@ class SourceInjectionConfig(Object):
optional=True)
f_max = Float.T(
default=1.,
help='Maximum wave frequency of interest to be used for source model '
'patch size scaling.',
optional=True)
@ -131,10 +142,13 @@ class SourceInjectionConfig(Object):
@property
def nucleation(self):
'''Nucleation points as 2D array
'''
Nucleation points as 2D array.
:returns: nucleation points to be injected
:rtype: :py:class:`numpy.ndarray`, shape: ``(N, 2)``
:returns:
Nucleation points to be injected.
:rtype:
:py:class:`numpy.ndarray`, shape: ``(N, 2)``
'''
if self.nucleation_x.shape != self.nucleation_y.shape:
raise IndexError(
@ -148,19 +162,30 @@ class SourceInjectionConfig(Object):
class SourceInjector(Object):
'''Extract sources to be injected in Grond
'''
Extract sources to be injected in Grond.
'''
config = SourceInjectionConfig.T(
default=SourceInjectionConfig())
def _extract_from_tsumaps(self, lat, lon):
'''
Get faults from tsumaps (NEAMTHM18) database for given location.
'''
limit = self.config.neamthm18_limit
kwargs = {} if limit is None else {limit.quantity: limit.target_value}
kwargs = {} if limit is None else {
limit.quantity: limit.get_target_value()}
fms, probs = tsumaps.get_probs(lat, lon, **kwargs)
logger.info('Extracted %d models from NEAMTHM18' % len(fms))
return list(fms[:, 0]), list(fms[:, 1]), list(fms[:, 2])
def _extract_from_faults(self, lat, lon, magnitude=None, mt=None):
'''
Get faults from chosen fault database for given location.
'''
conf = self.config
if conf.fault_radius_max is not None:
radius = conf.fault_radius_max
@ -224,6 +249,9 @@ class SourceInjector(Object):
lengths,
widths,
nparallel):
'''
Generate dynamic sources based on prior data and injection config.
'''
from pyrocko import parimap
@ -381,6 +409,9 @@ class SourceInjector(Object):
for src in sources for nx, ny in conf.nucleation]
def _generate_point_sources(self, event, strikes, dips, rakes):
'''
Generate point sources based on prior data and injection config.
'''
from pyrocko import moment_tensor as pmt
base_source = gf.MTSource.from_pyrocko_event(event)
@ -407,21 +438,41 @@ class SourceInjector(Object):
return sources
def get_sources(self, event, store, source_type='dynamic', nparallel=1):
'''Extract source models to be tested for the event location
:param event: Event to extract likely source models for. Both location,
magnitude and moment tensor information (if available) are used.
:type event: :py:class:`pyrocko.model.event.Event`
:param store: green's function store needed for source generation
:type store: :py:class:`pyrocko.gf.store.Store`
:param nparallel: Number of processes which run in parallel.
:type nparallel: int
:returns: List of source models to be injected.
:rtype: list of
:py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
'''
Extract source models to be tested for the event location.
:param event:
Event to extract likely source models for. Both location, magnitude
and moment tensor information (if available) are used.
:type event:
:py:class:`pyrocko.model.event.Event`
:param store:
Green's function store needed for source generation.
:type store:
:py:class:`pyrocko.gf.store.Store`
:param source_type:
If ``dynamic``, pseudo dynamic rupture models are extracted. If
``point``, moment tensor point sources are returned instead.
:type source_type:
optional, str
:param nparallel:
Number of processes which run in parallel.
:type nparallel:
optional, int
:returns:
Source models to be injected.
:rtype:
list of
:py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture` or
:py:class:`~pyrocko.gf.seismosizer.MTSource`
'''
logger.info('Extract sources based on injection configuration ...')
conf = self.config
lat, lon = event.lat, event.lon
@ -457,6 +508,9 @@ class SourceInjector(Object):
rakes %= 360.
rakes[rakes > 180.] -= 360.
logger.info('Generate %s sources for %d inputs ...' %
(source_type, len(strikes)))
if source_type == 'dynamic':
return self._generate_dynamic_sources(
event,
@ -481,7 +535,9 @@ def expand(base_path, path):
class Inverter(Object):
'''
Inverter base class.
'''
def run(self):
pass
@ -528,10 +584,13 @@ class GrondInverter(Inverter):
@property
def store(self):
'''Pyrocko Green's Function store
'''
Pyrocko Green's Function store.
:returns: green's function store
:rtype: :py:class:`pyrocko.gf.store.Store`
:returns:
Green's function store.
:rtype:
:py:class:`pyrocko.gf.store.Store`
'''
if self._store is not None:
return self._store
@ -548,10 +607,13 @@ class GrondInverter(Inverter):
@property
def event(self):
'''Pyrocko event
'''
Pyrocko event.
:returns: event
:rtype: :py:class:`pyrocko.model.event.Event`
:returns:
Event.
:rtype:
:py:class:`pyrocko.model.event.Event`
'''
if self._event is not None:
return self._event
@ -570,10 +632,10 @@ class GrondInverter(Inverter):
if isinstance(p, list):
for path in p:
if not op.exists(path):
logger.warn('Path %s does not exist.' % path)
logger.warning('Path %s does not exist.' % path)
else:
if not op.exists(p):
logger.warn('Path %s does not exist.' % p)
logger.warning('Path %s does not exist.' % p)
return p
@ -608,10 +670,13 @@ class GrondInverter(Inverter):
@property
def reportconf_path(self):
'''Path to the Grond report configuration file
'''
Path to the Grond report configuration file.
:returns: path to Grond report configuration file
:rtype: str
:returns:
Path to Grond report configuration file.
:rtype:
str
'''
if self.report_path is None:
raise ValueError('report_path needed')
@ -622,10 +687,13 @@ class GrondInverter(Inverter):
@property
def sources_inject(self):
'''Injected models from :py:class:`ewrica.si.invert.SourceInjector`
'''
Injected models from :py:class:`ewrica.si.invert.SourceInjector`.
:returns: sources to be injected within the run
:rtype: list of :py:class:`pyrocko.gf.seismosizer.Source`
:returns:
Sources to be injected within the run.
:rtype:
list of :py:class:`pyrocko.gf.seismosizer.Source`
'''
from grond.problems.cmt.problem import CMTProblemConfig
from grond.problems.dynamic_rupture.problem import \
@ -659,10 +727,13 @@ class GrondInverter(Inverter):
@property
def fn_sources_inject(self):
'''Path to the Grond source injection file
'''
Path to the Grond source injection file.
:returns: path to Grond source injection file
:rtype: str
:returns:
Path to Grond source injection file.
:rtype:
str
'''
if self.report_path is None:
raise ValueError('report_path needed')
@ -673,14 +744,16 @@ class GrondInverter(Inverter):
'sources_inject.yaml'])
def sources_to_injector_config(self):
'''Update Grond InjectionConfig for injected sources
'''
Update Grond InjectionConfig for injected sources.
'''
if self.grond_config is not None:
conf = self.grond_config.optimiser_config
sp = conf.sampler_phases
conf.sampler_phases = [InjectionSamplerPhase(
niterations=len(self.sources_inject),
sources_path=self.fn_sources_inject)]
sources_path=self.fn_sources_inject)] + sp
self.grond_config.optimiser_config = conf
@ -690,7 +763,8 @@ class GrondInverter(Inverter):
logger.warning('Could not update config - No config given')
def save_sources_inject(self):
'''Save inject sources in file
'''
Save inject sources in file.
'''
if self.fn_sources_inject is not None:
dump_all(self.sources_inject, filename=self.fn_sources_inject)
@ -700,7 +774,8 @@ class GrondInverter(Inverter):
logger.warning('Could not save sources - No filename given')
def generate_reportconf(self):
'''Generate and save Grond report config file
'''
Generate and save Grond report config file.
'''
if self.report_path is not None:
report_conf = report.ReportConfig(
@ -713,10 +788,13 @@ class GrondInverter(Inverter):
'Dumped specific grond report config %s', self.reportconf_path)
def run(self, force=False):
'''Run Grond optimization
'''
Run Grond optimization.
:param force: If True, possible existing older runs are overwritten
:type force: bool
:param force:
If ``True``, possible existing older runs are overwritten
:type force:
bool
'''
if self.grond_config is None:
raise ValueError('Config file needed to run Grond Inversion')
@ -731,7 +809,7 @@ class GrondInverter(Inverter):
os.chdir(self.base_path)
env = Environment(config=conf, event_names=[self.event_id])
env = Environment(config=conf, event_names=conf.event_names)
core.go(
environment=env,
@ -743,7 +821,8 @@ class GrondInverter(Inverter):
os.chdir(self._cwd)
def report(self):
'''Generate Grond report for run for easier result validation
'''
Generate Grond report for run for easier result validation.
'''
os.chdir(self.base_path)
@ -769,16 +848,18 @@ class GrondInverter(Inverter):
os.chdir(self._cwd)
def run_full(self, create_report=True, *args, **kwargs):
'''Run full optimization loop
'''
Run full optimization loop.
Different steps as source injection, report config generation, the
Grond run and reporting (if wished) are performed.
Arguments and keyword arguments for
:py:meth:`ewrica.si.invert.GrondInverter.run` can be given.
:param create_report: If True, a Grond report is created after
completing the run
:type create_report: optional, bool
:param create_report:
If ``True``, a Grond report is created after completing the run.
:type create_report:
optional, bool
'''
self.sources_to_injector_config()
self.save_sources_inject()

103
src/si/module.py

@ -2,12 +2,12 @@
#
# The Developers, 21st Century
import logging
# import os.path as op
import os.path as op
# from grond.environment import Environment
from ..module import BaseModule
from .invert import GrondInverter
# from .invert import GrondInverter
logger = logging.getLogger('ewrica.si.module')
@ -31,34 +31,93 @@ class SourceInversionModule(BaseModule):
# self.config.check_paths()
# TODO Check for event file, meta data and data, config files, paths
def initialize_inverter(self):
inv = {}
inv['grond'] = GrondInverter()
def run_point_source_inversion(
self,
problem,
force=False,
create_report=False):
self._inverter = inv
conf = self.config.pointsource_config.get_config()
def run_point_source_inversion(self, problem, force=False):
logger.info('running point source inversion')
if not conf.run:
logger.info('Point source inversion skipped')
return
if not self._inverter:
self.initialize_inverter()
logger.info('running point source inversion ...')
inv = self._inverter['grond']
inv.config = self.config.pointsource.get_config()
inv.config.set_basepath(problem.module_processing_dir(self))
store = problem.get_parent()
inv.event_id = problem.event_id
inv.run(force=force)
conf.grond_config_fn = op.join(
problem.config_dir, 'extra',
op.basename(conf.grond_config_fn))
conf.load_grond_config_from_file()
conf.grond_config.set_basepath('.', parent_path_prefix=store.store_dir)
conf.grond_config.path_prefix = store.store_dir
conf.grond_config.event_names = [problem.event_id]
conf.base_path = problem.module_processing_dir(self)
conf.report_path = problem.module_result_dir(self)
inv = conf.get_inverter(
problem_name='point_source',
event_id=problem.event_id)
dc = inv.grond_config.dataset_config
dc.set_basepath(
conf.grond_config.dataset_config.get_basepath(),
parent_path_prefix=store.store_dir)
inv.grond_config.dataset_config = dc
inv.run_full(force=force, create_report=create_report)
# TODO convert grond results to ewrica.sources
def run_ids_inversion(self):
logger.info('running IDS finite fault estimation')
pass
def run_dynamic_rupture_inversion(self):
logger.info('running pseudo dynamic ruputure inversion')
pass
def run_dynamic_rupture_inversion(
self,
problem,
force=False,
create_report=False):
conf = self.config.dynamic_rupture_config.get_config()
if not conf.run:
logger.info('Pseudo dynamic rupture inversion skipped')
return
logger.info('running pseudo dynamic ruputure inversion ...')
def run(self, problem=None, force=False):
store = problem.get_parent()
conf.grond_config_fn = op.join(
problem.config_dir, 'extra',
op.basename(conf.grond_config_fn))
conf.load_grond_config_from_file()
conf.grond_config.set_basepath('.', parent_path_prefix=store.store_dir)
conf.grond_config.path_prefix = store.store_dir
conf.grond_config.event_names = [problem.event_id]
conf.base_path = problem.module_processing_dir(self)
conf.report_path = problem.module_result_dir(self)
inv = conf.get_inverter(
problem_name='dynamic_rupture',
event_id=problem.event_id)
dc = inv.grond_config.dataset_config
dc.set_basepath(
conf.grond_config.dataset_config.get_basepath(),
parent_path_prefix=store.store_dir)
inv.grond_config.dataset_config = dc
inv.run_full(force=force, create_report=create_report)
# TODO convert grond results to ewrica.sources
def run(self, problem=None, *args, **kwargs):
# TODO call grond
# TODO use grond results (raw, so misfit and models) for statistics
# TODO convert grond results to own source models
@ -74,10 +133,10 @@ class SourceInversionModule(BaseModule):
logger.info('running source inversion module in "%s" mode' % mode)
self.run_point_source_inversion(problem, force=force)
self.run_point_source_inversion(problem, *args, **kwargs)