Compare commits

...

7 Commits

  1. 1
      setup.py
  2. 9
      src/apps/grond.py
  3. 28
      src/report/app/js/app.js
  4. 31
      src/scenario.py
  5. 1
      src/targets/__init__.py
  6. 6
      src/targets/base.py
  7. 1
      src/targets/phase_pick/__init__.py
  8. 128
      src/targets/phase_pick/plot.py
  9. 206
      src/targets/phase_pick/target.py
  10. 168
      test/test_locate.py

1
setup.py

@ -106,6 +106,7 @@ setup(
'grond.targets.waveform_oac',
'grond.targets.satellite',
'grond.targets.gnss_campaign',
'grond.targets.phase_pick',
'grond.problems',
'grond.problems.cmt',
'grond.problems.double_dc',

9
src/apps/grond.py

@ -418,6 +418,12 @@ def command_scenario(args):
'--gf-store-superdirs',
dest='gf_store_superdirs',
help='Comma-separated list of directories containing GF stores')
parser.add_option(
'--no-map',
dest='make_map',
default=True,
action='store_false',
help='suppress generation of map')
parser, options, args = cl_parse('scenario', args, setup)
@ -468,7 +474,8 @@ def command_scenario(args):
scenario.build(
force=options.force,
interactive=True,
gf_store_superdirs=gf_store_superdirs)
gf_store_superdirs=gf_store_superdirs,
make_map=options.make_map)
logger.info(CLIHints('scenario',
config=scenario.get_grond_config_path(),

28
src/report/app/js/app.js

@ -112,6 +112,8 @@ var yaml_type_map = [
['!grond.SatelliteTargetGroup', Dummy],
['!grond.PhaseRatioTarget', Dummy],
['!grond.PhaseRatioTargetGroup', Dummy],
['!grond.PhasePickTarget', Dummy],
['!grond.PhasePickTargetGroup', Dummy],
['!grond.FeatureMeasure', Dummy],
['!grond.CMTProblem', Dummy],
['!grond.RectangularProblem', Dummy],
@ -562,15 +564,17 @@ angular.module('reportApp', ['ngRoute', 'ngSanitize'])
YamlDoc.query(
get_path(problem_name) + '/stats.yaml',
function(doc) {
doc.title = 'Parameter Results';
doc.name = 'parameter results';
doc.variant = 'default';
doc.section = 'run';
doc.stats_path = get_path(problem_name) + '/stats.yaml';
doc.feather_icon = 'book';
doc.template = 'parameter-table';
insert_group(problem_name, doc);
if (doc) {
doc.title = 'Parameter Results';
doc.name = 'parameter results';
doc.variant = 'default';
doc.section = 'run';
doc.stats_path = get_path(problem_name) + '/stats.yaml';
doc.feather_icon = 'book';
doc.template = 'parameter-table';
insert_group(problem_name, doc);
}
},
{schema: report_schema});
@ -768,7 +772,11 @@ angular.module('reportApp', ['ngRoute', 'ngSanitize'])
};
$scope.doc_matches_keyword = function(doc) {
return ($filter('filter')(doc.items, $scope.keyword)).length > 0;
if (doc) {
return ($filter('filter')(doc.items, $scope.keyword)).length > 0;
} else {
return false;
}
};
})

31
src/scenario.py

@ -4,6 +4,7 @@ import shutil
import os.path as op
from pyrocko import gf, scenario, util
from pyrocko.gui import marker as pmarker
import math
import grond
@ -95,7 +96,12 @@ class GrondScenario(object):
for obs in self.observations],
source_generator=self.problem.get_scenario_source_generator())
def create_scenario(self, interactive=True, gf_store_superdirs=None):
def create_scenario(
self,
interactive=True,
gf_store_superdirs=None,
make_map=True):
logger.info('Creating scenario...')
scenario = self.get_scenario()
@ -128,11 +134,18 @@ class GrondScenario(object):
scenario.dump(filename=op.join(data_dir, 'scenario.yml'))
scenario.dump_data(path=data_dir)
scenario.make_map(op.join(self.project_dir, 'scenario_map.pdf'))
if make_map:
scenario.make_map(op.join(self.project_dir, 'scenario_map.pdf'))
shutil.move(op.join(data_dir, 'sources.yml'),
op.join(data_dir, 'scenario_sources.yml'))
markers = scenario.get_onsets()
marker_path = op.join(data_dir, 'picks', 'picks.markers')
if markers:
util.ensuredirs(marker_path)
pmarker.save_markers(markers, marker_path)
def get_grond_config(self):
engine_config = grond.EngineConfig(
gf_stores_from_pyrocko_config=False,
@ -163,13 +176,21 @@ class GrondScenario(object):
util.ensuredirs(config_path)
grond.write_config(self.get_grond_config(), config_path)
def build(self, force=False, interactive=False, gf_store_superdirs=None):
def build(
self,
force=False,
interactive=False,
gf_store_superdirs=None,
make_map=True):
logger.info('Building scenario...')
self.create_project_dir(force)
self.create_scenario(
interactive=interactive, gf_store_superdirs=gf_store_superdirs)
interactive=interactive,
gf_store_superdirs=gf_store_superdirs,
make_map=make_map)
self.create_grond_files()
@ -212,6 +233,8 @@ class WaveformObservation(Observation):
station_generator=scenario.targets.RandomStationGenerator(
nstations=self.nstations),
store_id=self.store_id,
tabulated_phases_from_store=True,
tabulated_phases_noise_scale=0.3,
seismogram_quantity='displacement')
def get_grond_target_group(self):

1
src/targets/__init__.py

@ -4,3 +4,4 @@ from .waveform_phase_ratio import * # noqa
from .waveform_oac import * # noqa
from .satellite import * # noqa
from .gnss_campaign import * # noqa
from .phase_pick import * # noqa

6
src/targets/base.py

@ -144,7 +144,11 @@ class MisfitTarget(Object):
def get_combined_weight(self):
if self._combined_weight is None:
self._combined_weight = num.ones(1, dtype=num.float)
w = self.manual_weight
for analyser in self.analyser_results.values():
w *= analyser.weight
self._combined_weight = num.array([w], dtype=num.float)
return self._combined_weight
def get_correlated_weights(self, nthreads=0):

1
src/targets/phase_pick/__init__.py

@ -0,0 +1 @@
from .target import * # noqa

128
src/targets/phase_pick/plot.py

@ -0,0 +1,128 @@
import numpy as num
from matplotlib import pyplot as plt, cm as mcm, colors as mcolors
from pyrocko.guts import Tuple, Float
from pyrocko.plot import mpl_init, mpl_margins
from grond.plot.config import PlotConfig
from grond.plot.collection import PlotItem
from .target import PhasePickTarget
class PickResidualsPlot(PlotConfig):
'''Travel time residuals plot.'''
name = 'pick_residuals'
size_cm = Tuple.T(
2, Float.T(),
default=(15., 7.5),
help='Width and length of the figure in [cm]')
def make(self, environ):
environ.setup_modelling()
cm = environ.get_plot_collection_manager()
mpl_init(fontsize=self.font_size)
ds = environ.get_dataset()
history = environ.get_history(subset='harvest')
optimiser = environ.get_optimiser()
cm.create_group_mpl(
self,
self.draw_figures(ds, history, optimiser),
title=u'Pick Residuals',
section='fits',
feather_icon='activity',
description=u'''
Travel time residuals.
''')
def draw_figures(self, ds, history, optimiser):
# gms = history.get_sorted_primary_misfits()[::-1]
models = history.get_sorted_primary_models()[::-1]
problem = history.problem
targets = [
t for t in problem.targets if
isinstance(t, PhasePickTarget)]
targets.sort(key=lambda t: t.distance_to(problem.base_source))
tpaths = sorted(set(t.path for t in targets))
ntargets = len(targets)
nmodels = history.nmodels
tts = num.zeros((nmodels, ntargets, 2))
# distances = num.zeros((nmodels, ntargets))
for imodel in range(nmodels):
model = models[imodel, :]
# source = problem.get_source(model)
results = problem.evaluate(model, targets=targets)
for itarget, result in enumerate(results):
result = results[itarget]
tts[imodel, itarget, :] = result.tobs, result.tsyn
# distances[imodel, itarget] = source.distance_to(
# targets[itarget])
ok = num.all(num.isfinite(tts), axis=2)
ok = num.all(ok, axis=0)
print(ok.size, ntargets)
iorder = num.arange(nmodels)
icolor = iorder
cmap = mcm.ScalarMappable(
norm=mcolors.Normalize(vmin=num.min(icolor), vmax=num.max(icolor)),
cmap=plt.get_cmap('coolwarm'))
imodel_to_color = []
for imodel in range(nmodels):
imodel_to_color.append(cmap.to_rgba(icolor[imodel]))
for tpath in tpaths:
mask = num.array(
[t.path == tpath for t in targets], dtype=num.bool)
mask &= ok
targets_this = [t for (it, t) in enumerate(targets) if mask[it]]
fontsize = self.font_size
item = PlotItem(
name='fig_%s' % tpath)
item.attributes['targets'] = [t.string_id() for t in targets_this]
print(item.attributes['targets'])
fig = plt.figure(figsize=self.size_inch)
labelpos = mpl_margins(
fig, nw=1, nh=1, left=7., right=1., bottom=10., top=3,
units=fontsize)
axes = fig.add_subplot(1, 1, 1)
labelpos(axes, 2.5, 2.0)
for imodel in range(0, nmodels, 10):
color = imodel_to_color[imodel]
axes.plot(
num.arange(len(targets_this)),
# distances[imodel, mask],
# tts[imodel, mask, 0],
tts[imodel, mask, 1] - tts[imodel, mask, 0],
'_',
ms=2.0,
color=color)
yield item, fig
def get_plot_classes():
return [
PickResidualsPlot]

206
src/targets/phase_pick/target.py

@ -0,0 +1,206 @@
import logging
import numpy as num
from pyrocko.guts import Float, Tuple, String
from pyrocko import gf
from ..base import (
MisfitTarget, TargetGroup, MisfitResult)
from grond.meta import has_get_plot_classes
guts_prefix = 'grond'
logger = logging.getLogger('grond.targets.phase_pick.target')
def log_exclude(target, reason):
logger.debug('Excluding potential target %s: %s' % (
target.string_id(), reason))
class PhasePickTargetGroup(TargetGroup):
'''
Generate targets to fit phase arrival times.
'''
distance_min = Float.T(optional=True)
distance_max = Float.T(optional=True)
distance_3d_min = Float.T(optional=True)
distance_3d_max = Float.T(optional=True)
depth_min = Float.T(optional=True)
depth_max = Float.T(optional=True)
store_id = gf.StringID.T(optional=True)
pick_synthetic_traveltime = gf.Timing.T(
help='Synthetic phase arrival definition.')
pick_phasename = String.T(
help='Name of phase in pick file.')
def get_targets(self, ds, event, default_path='none'):
logger.debug('Selecting phase pick targets...')
origin = event
targets = []
for st in ds.get_stations():
target = PhasePickTarget(
codes=st.nsl(),
lat=st.lat,
lon=st.lon,
north_shift=st.north_shift,
east_shift=st.east_shift,
depth=st.depth,
store_id=self.store_id,
manual_weight=self.weight,
normalisation_family=self.normalisation_family,
path=self.path or default_path,
pick_synthetic_traveltime=self.pick_synthetic_traveltime,
pick_phasename=self.pick_phasename)
if self.distance_min is not None and \
target.distance_to(origin) < self.distance_min:
log_exclude(target, 'distance < distance_min')
continue
if self.distance_max is not None and \
target.distance_to(origin) > self.distance_max:
log_exclude(target, 'distance > distance_max')
continue
if self.distance_3d_min is not None and \
target.distance_3d_to(origin) < self.distance_3d_min:
log_exclude(target, 'distance_3d < distance_3d_min')
continue
if self.distance_3d_max is not None and \
target.distance_3d_to(origin) > self.distance_3d_max:
log_exclude(target, 'distance_3d > distance_3d_max')
continue
if self.depth_min is not None and \
target.depth < self.depth_min:
log_exclude(target, 'depth < depth_min')
continue
if self.depth_max is not None and \
target.depth > self.depth_max:
log_exclude(target, 'depth > depth_max')
continue
marker = ds.get_pick(
event.name,
target.codes[:3],
target.pick_phasename)
if not marker:
log_exclude(target, 'no pick available')
continue
target.set_dataset(ds)
targets.append(target)
return targets
class PhasePickResult(MisfitResult):
tobs = Float.T(optional=True)
tsyn = Float.T(optional=True)
@has_get_plot_classes
class PhasePickTarget(gf.Location, MisfitTarget):
'''
Target to fit phase arrival times.
'''
codes = Tuple.T(
3, String.T(),
help='network, station, location codes.')
store_id = gf.StringID.T(
help='ID of Green\'s function store (only used for earth model).')
pick_synthetic_traveltime = gf.Timing.T(
help='Synthetic phase arrival definition.')
pick_phasename = String.T(
help='Name of phase in pick file.')
can_bootstrap_weights = True
def __init__(self, **kwargs):
gf.Location.__init__(self, **kwargs)
MisfitTarget.__init__(self, **kwargs)
self._tobs_cache = {}
@classmethod
def get_plot_classes(cls):
from . import plot
plots = super(PhasePickTarget, cls).get_plot_classes()
plots.extend(plot.get_plot_classes())
return plots
def string_id(self):
return '.'.join(x for x in (self.path,) + self.codes)
def set_dataset(self, ds):
MisfitTarget.set_dataset(self, ds)
self._tobs_cache = {}
def get_plain_targets(self, engine, source):
return self.prepare_modelling(engine, source, None)
def prepare_modelling(self, engine, source, targets):
return []
def get_times(self, engine, source):
tsyn = None
k = source.name
if k not in self._tobs_cache:
ds = self.get_dataset()
tobs = None
marker = ds.get_pick(
source.name,
self.codes[:3],
self.pick_phasename)
if marker:
self._tobs_cache[k] = marker.tmin
else:
self._tobs_cache[k] = None
tobs = self._tobs_cache[k]
store = engine.get_store(self.store_id)
tsyn = source.time + store.t(
self.pick_synthetic_traveltime, source, self)
return tobs, tsyn
def finalize_modelling(
self, engine, source, modelling_targets, modelling_results):
ds = self.get_dataset() # noqa
tobs, tsyn = self.get_times(engine, source)
if None not in (tobs, tsyn):
misfit = abs(tobs - tsyn)
else:
misfit = None
norm = 1.0
result = PhasePickResult(
tobs=tobs,
tsyn=tsyn,
misfits=num.array([[misfit, norm]], dtype=num.float))
return result
__all__ = '''
PhasePickTargetGroup
PhasePickTarget
PhasePickResult
'''.split()

168
test/test_locate.py

@ -0,0 +1,168 @@
from __future__ import absolute_import
import os
import shutil
from pyrocko import guts
from grond import config
from grond import Environment
from . import common
from .common import grond, chdir
_multiprocess_can_split = True
def test_locate():
playground_dir = common.get_playground_dir()
common.get_test_data('gf_stores/crust2_ib/')
gf_stores_path = common.test_data_path('gf_stores')
with chdir(playground_dir):
scenario_dir = 'scenario_locate'
if os.path.exists(scenario_dir):
shutil.rmtree(scenario_dir)
grond('scenario', '--targets=waveforms', '--nevents=1',
'--nstations=6', '--gf-store-superdirs=%s' % gf_stores_path,
'--no-map',
scenario_dir)
with chdir(scenario_dir):
config_path = 'config/scenario.gronf'
quick_config_path = 'config/scenario_quick.gronf'
event_names = grond('events', config_path).strip().split('\n')
env = Environment([config_path] + event_names)
conf = env.get_config()
mod_conf = conf.clone()
target_groups = guts.load_all(string='''
--- !grond.PhasePickTargetGroup
normalisation_family: picks
path: pick.p
weight: 1.0
store_id: crust2_ib
distance_min: 10000.0
distance_max: 1000000.0
pick_synthetic_traveltime: '{stored:any_P}'
pick_phasename: 'any_P'
--- !grond.PhasePickTargetGroup
normalisation_family: picks
path: pick.s
weight: 1.0
store_id: crust2_ib
distance_min: 10000.0
distance_max: 1000000.0
pick_synthetic_traveltime: '{stored:any_S}'
pick_phasename: 'any_S'
''')
mod_conf.dataset_config.picks_paths = [
'data/scenario/picks/picks.markers']
mod_conf.target_groups = target_groups
# mod_conf.set_elements(
# 'analyser_configs[:].niterations', 100)
# mod_conf.set_elements(
# 'optimiser_config.sampler_phases[:].niterations', 100)
# mod_conf.set_elements(
# 'optimiser_config.nbootstrap', 5)
mod_conf.optimiser_config.sampler_phases[-1].niterations = 3000
mod_conf.set_basepath(conf.get_basepath())
config.write_config(mod_conf, quick_config_path)
grond('diff', config_path, quick_config_path)
grond('check', quick_config_path, *event_names)
grond('go', quick_config_path, *event_names)
# rundir_paths = common.get_rundir_paths(
# quick_config_path, event_names)
# grond('report', *rundir_paths)
def test_joint_locate():
playground_dir = common.get_playground_dir()
common.get_test_data('gf_stores/crust2_ib/')
gf_stores_path = common.test_data_path('gf_stores')
with chdir(playground_dir):
scenario_dir = 'scenario_joint_locate'
if os.path.exists(scenario_dir):
shutil.rmtree(scenario_dir)
grond('scenario', '--targets=waveforms', '--nevents=1',
'--nstations=6', '--gf-store-superdirs=%s' % gf_stores_path,
'--no-map',
scenario_dir)
with chdir(scenario_dir):
config_path = 'config/scenario.gronf'
quick_config_path_templ = 'config/scenario_quick_%i.gronf'
event_names = grond('events', config_path).strip().split('\n')
env = Environment([config_path] + event_names)
conf = env.get_config()
for irun, pick_weight in enumerate([0.5, 1.0, 2.0]):
quick_config_path = quick_config_path_templ % irun
mod_conf = conf.clone()
target_groups = guts.load_all(string='''
--- !grond.PhasePickTargetGroup
normalisation_family: picks
path: pick.p
weight: 1.0
store_id: crust2_ib
distance_min: 10000.0
distance_max: 1000000.0
pick_synthetic_traveltime: '{stored:any_P}'
pick_phasename: 'any_P'
--- !grond.PhasePickTargetGroup
normalisation_family: picks
path: pick.s
weight: 1.0
store_id: crust2_ib
distance_min: 10000.0
distance_max: 1000000.0
pick_synthetic_traveltime: '{stored:any_S}'
pick_phasename: 'any_S'
''')
mod_conf.problem_config.name_template \
= 'cmt_${event_name}_%i' % irun
mod_conf.dataset_config.picks_paths = [
'data/scenario/picks/picks.markers']
mod_conf.target_groups.extend(target_groups)
mod_conf.set_elements(
'target_groups[-2:].weight', pick_weight)
# mod_conf.set_elements(
# 'analyser_configs[:].niterations', 100)
# mod_conf.set_elements(
# 'optimiser_config.sampler_phases[:].niterations', 100)
# mod_conf.set_elements(
# 'optimiser_config.nbootstrap', 5)
mod_conf.optimiser_config.sampler_phases[-1].niterations \
= 10000
mod_conf.set_elements(
'optimiser_config.sampler_phases[:].seed', 23)
mod_conf.set_basepath(conf.get_basepath())
config.write_config(mod_conf, quick_config_path)
grond('diff', config_path, quick_config_path)
grond('check', quick_config_path, *event_names)
grond('go', quick_config_path, *event_names)
rundir_paths = common.get_rundir_paths(
quick_config_path, event_names)
grond('report', *rundir_paths)
Loading…
Cancel
Save