Compare commits

...

11 Commits

Author SHA1 Message Date
Pyrocko Tester 1f43848b54 pick residual plots 3 months ago
Pyrocko Tester fdf9caa2ad grond export: add --effective-lat-lon option 3 months ago
Sebastian Heimann 85b57fcfbe locate with picks: pick residuals plot (wip) 1 year ago
Sebastian Heimann d0bb833f6d use String instead of StringID for event names 3 months ago
Pyrocko Tester 60c5aa4897 report: respect filter when clicking "Select all" 1 year ago
Sebastian Heimann 93bd5782ef locating with picks: fix weighting issue 1 year ago
Pyrocko Tester 3e66ff26be picks: handle in report 1 year ago
Pyrocko Tester 514bb30a88 picks: improve robustness 1 year ago
Pyrocko Tester a69e2668c5 locate with picks: handle incomplete picksets, cache observed picks 2 years ago
Sebastian Heimann 5018bd9e7e testing: improve speed of locate test 2 years ago
Sebastian Heimann 659eb42309 pick residual fitting (Grond can now be used as a locator) 2 years ago
  1. 1
      setup.py
  2. 18
      src/apps/grond.py
  3. 6
      src/config.py
  4. 10
      src/core.py
  5. 5
      src/plot/config.py
  6. 46
      src/report/app/js/app.js
  7. 2
      src/report/app/templates/report_list_modal.tmpl.html
  8. 31
      src/scenario.py
  9. 1
      src/targets/__init__.py
  10. 6
      src/targets/base.py
  11. 12
      src/targets/gnss_campaign/plot.py
  12. 1
      src/targets/phase_pick/__init__.py
  13. 373
      src/targets/phase_pick/plot.py
  14. 206
      src/targets/phase_pick/target.py
  15. 131
      src/targets/plot.py
  16. 12
      src/targets/waveform/plot.py
  17. 168
      test/test_locate.py

1
setup.py

@ -134,6 +134,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',

18
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(),
@ -1016,6 +1023,12 @@ def command_export(args):
'--output', dest='filename', metavar='FILE',
help='write output to FILE')
parser.add_option(
'--effective-lat-lon', dest='effective_lat_lon',
action='store_true',
help='convert north_shift/east_shift offsets to true lat/lon '
'coordinates (when outputting event objects).')
parser, options, args = cl_parse('export', args, setup)
if len(args) < 2:
help_and_die(parser, 'arguments required')
@ -1044,7 +1057,8 @@ def command_export(args):
filename=options.filename,
type=options.type,
pnames=pnames,
selection=options.selection)
selection=options.selection,
effective_lat_lon=options.effective_lat_lon)
except grond.GrondError as e:
die(str(e))

6
src/config.py

@ -1,7 +1,7 @@
import re
import os.path as op
from pyrocko import guts, gf
from pyrocko.guts import Bool, List
from pyrocko.guts import Bool, List, String
from .meta import Path, HasPaths, GrondError
from .dataset import DatasetConfig
@ -82,11 +82,11 @@ class Config(HasPaths):
default=EngineConfig.D(),
help=':class:`pyrocko.gf.LocalEngine` configuration')
event_names = List.T(
gf.StringID.T(),
String.T(),
help='Restrict application to given event names. If empty, all events '
'found through the dataset configuration are considered.')
event_names_exclude = List.T(
gf.StringID.T(),
String.T(),
help='Event names to be excluded')
def __init__(self, *args, **kwargs):

10
src/core.py

@ -688,7 +688,8 @@ def format_stats(rs, fmt):
def export(
what, rundirs, type=None, pnames=None, filename=None, selection=None):
what, rundirs, type=None, pnames=None, filename=None, selection=None,
effective_lat_lon=False):
if pnames is not None:
pnames_clean = [pname.split('.')[0] for pname in pnames]
@ -729,14 +730,21 @@ def export(
elif type == 'source':
source = problem.get_source(x)
if effective_lat_lon:
source.set_origin(*source.effective_latlon)
guts.dump(source, stream=out)
elif type == 'event':
ev = problem.get_source(x).pyrocko_event()
if effective_lat_lon:
ev.set_origin(*ev.effective_latlon)
model.dump_events([ev], stream=out)
elif type == 'event-yaml':
ev = problem.get_source(x).pyrocko_event()
if effective_lat_lon:
ev.set_origin(*ev.effective_latlon)
guts.dump_all([ev], stream=out)
else:

5
src/plot/config.py

@ -6,6 +6,7 @@ guts_prefix = 'grond'
inch = 2.54
points = inch / 72.0
class PlotFormat(Object):
@ -145,6 +146,10 @@ class PlotConfig(Object):
def size_inch(self):
return self.size_cm[0]/inch, self.size_cm[1]/inch
@property
def size_points(self):
return self.size_cm[0]/points, self.size_cm[1]/points
def make(self, environ):
pass

46
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.DoubleDCProblem', Dummy],
@ -291,6 +293,15 @@ angular.module('reportApp', ['ngRoute', 'ngSanitize'])
return selected_problem_names;
};
funcs.set_selected = function(problem_names) {
selected_problem_names.length = 0;
for (var i=0; i<problem_names.length; i++) {
selected_problem_names.push(problem_names[i]);
}
selected_problem_names.sort()
};
funcs.is_selected = function(problem_name) {
return includes(selected_problem_names, problem_name);
};
@ -420,6 +431,15 @@ angular.module('reportApp', ['ngRoute', 'ngSanitize'])
return $filter('filter')(ordered_lines[order_skey], $scope.list_search_keyword);
};
$scope.select_all_filtered = function() {
var lines = $scope.get_ordered_report_entries();
var problem_names = [];
for (var i=0; i<lines.length; i++) {
problem_names.push(
lines[i].report_entry.problem_name);
}
rl.set_selected(problem_names);
};
})
.controller('ReportController', function(
@ -564,15 +584,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});
@ -770,7 +792,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;
}
};
})

2
src/report/app/templates/report_list_modal.tmpl.html

@ -83,7 +83,7 @@
<button type="button" class="btn btn-secondary mr-2" ng-click="rl.reload()">
<span style="width: 1em; height: 1em; margin: 0.2em;" data-feather="refresh-cw"></span> Reload
</button>
<button type="button" class="btn btn-secondary mr-2" ng-click="rl.selected_all()">
<button type="button" class="btn btn-secondary mr-2" ng-click="select_all_filtered()">
Select all
</button>
<button type="button" class="btn btn-secondary mr-2" ng-click="rl.selected_none()">

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):

12
src/targets/gnss_campaign/plot.py

@ -277,22 +277,22 @@ components).
item = PlotItem(name='station_distribution-N-%s' % target.path)
fig, ax, legend = self.plot_station_distribution(
azimuths, distances, ws_n, labels)
legend.set_title('Weight, N components')
azimuths, distances, ws_n, labels,
legend_title='Weight, N components')
yield (item, fig)
item = PlotItem(name='station_distribution-E-%s' % target.path)
fig, ax, legend = self.plot_station_distribution(
azimuths, distances, ws_e, labels)
legend.set_title('Weight, E components')
azimuths, distances, ws_e, labels,
legend_title='Weight, E components')
yield (item, fig)
item = PlotItem(name='station_distribution-U-%s' % target.path)
fig, ax, legend = self.plot_station_distribution(
azimuths, distances, ws_u, labels)
legend.set_title('Weight, U components')
azimuths, distances, ws_u, labels,
legend_title='Weight, U components')
yield (item, fig)

1
src/targets/phase_pick/__init__.py

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

373
src/targets/phase_pick/plot.py

@ -0,0 +1,373 @@
import math
from matplotlib import pyplot as plt, colors as mcolors
import numpy as num
from pyrocko import util
from pyrocko.guts import Tuple, Float
from pyrocko.plot import mpl_init, mpl_margins
from grond import meta
from grond.plot.config import PlotConfig
from grond.plot.collection import PlotItem
from .target import PhasePickTarget
from ..plot import StationDistributionPlot
km = 1000.
d2r = math.pi / 180.
class PickResidualsPlot(PlotConfig):
'''Travel time residuals plot.'''
name = 'pick_residuals'
size_cm = Tuple.T(
2, Float.T(),
default=(15.0, 10.0),
help='Width and height 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='watch',
description=u'''
Travel time residuals.
The difference between observed and synthetic phase arrival times is shown
for a subset of the models in the solution ensemble. Models with good global
performance are shown in red, less good ones in blue. Vertical black bars mark
the classical best model.
''')
def draw_figures(self, ds, history, optimiser):
show_mean_residuals = False
# 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))
azimuths = 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])
azimuths[imodel, itarget] = source.azibazi_to(
targets[itarget])[0]
ok = num.all(num.isfinite(tts), axis=2)
ok = num.all(ok, axis=0)
targets = [
target for (itarget, target) in enumerate(targets) if ok[itarget]]
tts = tts[:, ok, :]
distances = distances[:, ok]
azimuths = azimuths[:, ok]
residuals = tts[:, :, 0] - tts[:, :, 1]
residual_amax = num.max(num.abs(residuals))
mean_residuals = num.mean(residuals, axis=0)
mean_distances = num.mean(distances, axis=0)
isort = num.array(
[x[2] for x in sorted(
zip(targets, mean_distances, range(len(targets))),
key=lambda x: (x[0].path, x[1]))], dtype=num.int)
distances = distances[:, isort]
distance_min = num.min(distances)
distance_max = num.max(distances)
azimuths = azimuths[:, isort]
targets = [targets[i] for i in isort]
ntargets = len(targets)
residuals = residuals[:, isort]
mean_residuals = mean_residuals[isort]
icolor = num.arange(nmodels)
norm = mcolors.Normalize(vmin=num.min(icolor), vmax=num.max(icolor))
cmap = plt.get_cmap('coolwarm')
napprox = max(
1, int(round((self.size_points[1] / self.font_size - 7))))
npages = (ntargets-1) // napprox + 1
ntargets_page = (ntargets-1) // npages + 1
for ipage in range(npages):
ilo, ihi = ipage*ntargets_page, (ipage+1)*ntargets_page
ihi = min(len(targets), ihi)
targets_page = targets[ilo:ihi]
item = PlotItem(
name='fig_%02i' % ipage)
item.attributes['targets'] = [t.string_id() for t in targets_page]
fig = plt.figure(figsize=self.size_inch)
labelpos = mpl_margins(
fig, nw=2, nh=1, left=12., right=1., bottom=5., top=1.,
wspace=1., units=self.font_size)
axes1 = fig.add_subplot(1, 3, 1)
axes2 = fig.add_subplot(1, 3, 2)
axes3 = fig.add_subplot(1, 3, 3)
for axes in (axes1, axes2, axes3):
axes.set_ylim(-0.5, len(targets_page)-0.5)
axes.invert_yaxis()
labelpos(axes, 2.5, 2.0)
axes1.axvline(0.0, color='black')
labels = []
lastpath = None
for ipos, t in enumerate(targets_page):
scodes = '.'.join(t.codes)
if lastpath is None or lastpath != t.path:
labels.append(t.path + '.' + scodes)
if lastpath is not None:
for axes in (axes1, axes2, axes3):
axes.axhline(ipos-0.5, color='black')
else:
labels.append(scodes)
lastpath = t.path
for ipos, itarget in enumerate(range(ilo, ihi)):
axes1.plot(
residuals[-1, itarget],
ipos,
'|',
ms=self.font_size,
color='black')
if show_mean_residuals:
axes1.plot(
mean_residuals[itarget],
ipos,
's',
ms=self.font_size,
color='none',
mec='black')
axes1.scatter(
residuals[::10, itarget],
util.num_full(icolor[::10].size, ipos, dtype=num.float),
c=icolor[::10],
cmap=cmap,
norm=norm,
alpha=0.5)
axes2.scatter(
distances[::10, itarget]/km,
util.num_full(icolor[::10].size, ipos, dtype=num.float),
c=icolor[::10],
cmap=cmap,
norm=norm,
alpha=0.5)
axes3.scatter(
azimuths[::10, itarget],
util.num_full(icolor[::10].size, ipos, dtype=num.float),
c=icolor[::10],
cmap=cmap,
norm=norm,
alpha=0.5)
axes1.set_yticks(num.arange(len(labels)))
axes1.set_yticklabels(labels)
axes1.set_xlabel('$T_{obs} - T_{syn}$ [s]')
axes1.set_xlim(-residual_amax, residual_amax)
axes2.set_xlabel('Distance [km]')
axes2.set_xlim(distance_min/km, distance_max/km)
axes3.set_xlabel('Azimuth [deg]')
axes3.set_xlim(-180., 180.)
axes2.get_yaxis().set_ticks([])
axes2.get_yaxis().set_ticks([])
axes3.get_yaxis().set_ticks([])
axes3.get_yaxis().set_ticks([])
yield item, fig
class PickResidualsStationDistribution(StationDistributionPlot):
'''
Plot showing residuals of seismic phase arrivals by azimuth and distance.
'''
name = 'pick_residuals_stations'
def make(self, environ):
environ.setup_modelling()
cm = environ.get_plot_collection_manager()
mpl_init(fontsize=self.font_size)
problem = environ.get_problem()
dataset = environ.get_dataset()
history = environ.get_history(subset='harvest')
cm.create_group_mpl(
self,
self.draw_figures(problem, dataset, history),
title=u'Pick Residuals by Location',
section='fits',
feather_icon='target',
description=u'''
Plot showing residuals for seismic phase arrivals by azimuth and distance.
Station locations in dependence of distance and azimuth are shown. The center
of the plot corresponds to the origin of the search space. Optimized source
locations are shown as small black dots (subset of the solution ensemble).
''')
def draw_figures(self, problem, dataset, history):
target_index = {}
i = 0
for target in problem.targets:
target_index[target] = i, i+target.nmisfits
i += target.nmisfits
misfits = history.misfits[history.get_sorted_misfits_idx(), ...]
gcms = problem.combine_misfits(
misfits[:1, :, :], get_contributions=True)[0, :]
models = history.get_sorted_primary_models()[::-1]
origin = problem.base_source
targets = [
t for t in problem.targets if isinstance(t, PhasePickTarget)]
ntargets = len(targets)
nmodels = history.nmodels
tts = num.zeros((nmodels, ntargets, 2))
sdata = []
for imodel in range(nmodels):
model = models[imodel, :]
source = problem.get_source(model)
sdata.append(
(origin.distance_to(source), origin.azibazi_to(source)[0]))
results = problem.evaluate(model, targets=targets)
for itarget, result in enumerate(results):
result = results[itarget]
tts[imodel, itarget, :] = result.tobs, result.tsyn
ok = num.all(num.isfinite(tts), axis=2)
ok = num.all(ok, axis=0)
targets_ok = [
target for (itarget, target) in enumerate(targets) if ok[itarget]]
tts = tts[:, ok, :]
residuals = tts[:, :, 0] - tts[:, :, 1]
mean_residuals = num.mean(residuals, axis=0)
rlo, rhi = num.percentile(mean_residuals, [10., 90.])
residual_amax = max(abs(rlo), abs(rhi))
target_to_residual = dict(
(target, residual)
for (target, residual)
in zip(targets_ok, mean_residuals))
cg_to_targets = meta.gather(targets, lambda t: (t.path, ))
cgs = sorted(cg_to_targets.keys())
for cg in cgs:
cg_str = '.'.join(cg)
targets = cg_to_targets[cg]
if len(targets) == 0:
continue
assert all(target_index[target][0] == target_index[target][1] - 1
for target in targets)
itargets = num.array(
[target_index[target][0] for target in targets])
labels = ['.'.join(x for x in t.codes[:3] if x) for t in targets]
azimuths = num.array(
[origin.azibazi_to(t)[0] for t in targets])
distances = num.array(
[t.distance_to(origin) for t in targets])
residuals = num.array(
[target_to_residual.get(t, num.nan) for t in targets])
item = PlotItem(
name='picks_contributions_%s' % cg_str,
title=u'Pick residuals and contributions (%s)' % cg_str,
description=u'\n\nMarkers are scaled according to their '
u'misfit contribution for the globally best '
u'source model.')
fig, axes, legend = self.plot_station_distribution(
azimuths, distances,
gcms[itargets],
labels,
colors=residuals,
cnorm=(-residual_amax, residual_amax),
cmap='RdYlBu',
clabel='$T_{obs} - T_{syn}$ [s]',
scatter_kwargs=dict(alpha=1.0),
legend_title='Contribution')
sources = [
problem.get_source(x) for x in models[::10]]
azimuths = num.array(
[origin.azibazi_to(s)[0] for s in sources])
distances = num.array(
[s.distance_to(origin) for s in sources])
axes.plot(
azimuths*d2r, distances, 'o', ms=2.0, color='black', alpha=0.3)
yield (item, fig)
def get_plot_classes():
return [
PickResidualsPlot,
PickResidualsStationDistribution]

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()

131
src/targets/plot.py

@ -5,7 +5,7 @@ import math
from pyrocko import plot
from pyrocko.guts import Tuple, Float
from matplotlib import pyplot as plt
from matplotlib import pyplot as plt, colors as mcolors, cm as mcm
from matplotlib import lines
from matplotlib.ticker import FuncFormatter
@ -18,6 +18,33 @@ d2r = math.pi / 180.
logger = logging.getLogger('targets.plot')
def ndig_inc(vinc):
assert vinc != 0.0
ndig = -math.floor(math.log10(abs(vinc)))
if ndig <= 0:
return 0
else:
return ndig + len(('%5.3f' % (vinc * 10**ndig)).rstrip('0')) - 2
def make_scale(vmin, vmax, approx_ticks=3, mode='0-max', snap=True):
cscale = plot.AutoScaler(approx_ticks=approx_ticks, mode=mode, snap=snap)
vmin, vmax, vinc = cscale.make_scale((vmin, vmax))
imin = math.ceil(vmin/vinc)
imax = math.floor(vmax/vinc)
vs = num.arange(imin, imax+1) * vinc
vexp = cscale.make_exp(vinc)
vexp_factor = 10**vexp
vinc_base = vinc / vexp_factor
ndig = ndig_inc(vinc_base)
return vs, vexp, '%%.%if' % ndig
def darken(f, c):
c = num.asarray(c)
return f*c
class StationDistributionPlot(PlotConfig):
''' Plot showing all waveform fits for the ensemble of solutions'''
@ -35,7 +62,9 @@ class StationDistributionPlot(PlotConfig):
def plot_station_distribution(
self, azimuths, distances, weights, labels=None,
scatter_kwargs=dict(), annotate_kwargs=dict(), maxsize=10**2):
colors=None, cmap=None, cnorm=None, clabel=None,
scatter_kwargs=dict(), annotate_kwargs=dict(), maxsize=10**2,
legend_title=''):
invalid_color = plot.mpl_color('aluminium3')
@ -80,15 +109,42 @@ class StationDistributionPlot(PlotConfig):
if weights_ref == 0.:
weights_ref = 1.0
colors = [scatter_default['c'] if s else invalid_color
for s in valid]
if colors is None:
scolors = [
scatter_default['c'] if s else invalid_color for s in valid]
else:
if cnorm is None:
cnorm = mcolors.Normalize()
elif isinstance(cnorm, tuple):
cnorm = mcolors.Normalize(vmin=cnorm[0], vmax=cnorm[1])
if cmap is None:
cmap = plt.get_cmap('viridis')
elif isinstance(cmap, str):
cmap = plt.get_cmap(cmap)
sm = mcm.ScalarMappable(norm=cnorm, cmap=cmap)
sm.set_array(colors)
scolors = [
sm.to_rgba(value) for value in colors]
scolors = num.array(scolors)
scatter_default.pop('c')
weights_scaled = (weights / weights_ref) * maxsize
ws, exp, fmt = make_scale(0., weights_ref)
ws = ws[1:]
weight_clip_min = ws[0]
weight_clip_min_scaled = (weight_clip_min / weights_ref) * maxsize
weights_scaled = num.maximum(weight_clip_min_scaled, weights_scaled)
stations = ax.scatter(
azimuths*d2r, distances, s=weights_scaled, c=colors,
azimuths*d2r, distances, s=weights_scaled, c=scolors,
edgecolors=darken(0.5, scolors),
linewidths=1.0,
**scatter_default)
if len(labels) < 30: # TODO: remove after impl. of collision detection
@ -101,53 +157,56 @@ class StationDistributionPlot(PlotConfig):
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.tick_params('y', labelsize=self.font_size, labelcolor='gray')
ax.grid(alpha=.3)
ax.grid(alpha=.2)
ax.set_ylim(0, distances.max()*1.1)
ax.yaxis.set_major_locator(plt.MaxNLocator(4))
ax.yaxis.set_major_formatter(
FuncFormatter(lambda x, pos: '%d km' % (x/km)))
FuncFormatter(lambda x, pos: '%d km' % (x/km) if x != 0.0 else ''))
# Legend
entries = 4
valid_marker = num.argmax(valid)
ecl = stations.get_edgecolor()
fc = tuple(stations.get_facecolor()[valid_marker])
ec = tuple(ecl[min(valid_marker, len(ecl)-1)])
def get_min_precision(values):
sig_prec = num.floor(
num.isfinite(num.log10(values[values > 0])))
if sig_prec.size == 0:
return 1
return int(abs(sig_prec.min())) + 1
legend_artists = [
lines.Line2D(
[0], [0], ls='none',
marker='o', ms=num.sqrt(rad), mfc=fc, mec=ec)
for rad in num.linspace(maxsize, .1*maxsize, entries)
]
sig_prec = get_min_precision(weights)
legend_annot = [
'{value:.{prec}f}'.format(value=val, prec=sig_prec)
for val in num.linspace(weights_ref, .1*weights_ref, entries)
]
legend_data = []
for w in ws[::-1]:
legend_data.append((
lines.Line2D(
[0], [0],
ls='none',
marker='o',
ms=num.sqrt(w/weights_ref*maxsize),
mfc=fc,
mec=ec),
fmt % (w/10**exp) + (
'$\\times 10^{%i}$' % exp if exp != 0 else '')))
if not num.all(valid):
legend_artists.append(
legend_data.append((
lines.Line2D(
[0], [0], ls='none',
marker='o', ms=num.sqrt(maxsize),
mfc=invalid_color, mec=invalid_color))
legend_annot.append('Excluded')
[0], [0],
marker='o',
ms=num.sqrt(maxsize),
mfc=invalid_color,
mec=darken(0.5, invalid_color),
mew=1.0),
'Excluded'))
legend = fig.legend(
legend_artists, legend_annot,
fontsize=self.font_size, loc=4,
*zip(*legend_data),
fontsize=self.font_size,
loc='upper left', bbox_to_anchor=(0.77, 0.4),
markerscale=1, numpoints=1,
frameon=False)
legend.set_title(
legend_title,
prop=dict(size=self.font_size))
cb_axes = fig.add_axes([0.8, 0.6, 0.02, 0.3])
if clabel is not None:
fig.colorbar(sm, cax=cb_axes, label=clabel, extend='both')
return fig, ax, legend

12
src/targets/waveform/plot.py

@ -1319,10 +1319,8 @@ location of the source.
u'weighting factor of the corresponding target\'s '
u'contribution in the misfit function.')
fig, ax, legend = self.plot_station_distribution(
azimuths, distances, ws[itargets], labels)
legend.set_title(
'Weight',
prop=dict(size=self.font_size))
azimuths, distances, ws[itargets], labels,
legend_title='Weight')
yield (item, fig)
@ -1334,10 +1332,8 @@ location of the source.
u'misfit contribution for the globally best '
u'source model.')
fig, ax, legend = self.plot_station_distribution(
azimuths, distances, gcms[itargets], labels)
legend.set_title(
'Contribution',
prop=dict(size=self.font_size))
azimuths, distances, gcms[itargets], labels,
legend_title='Contribution')
yield (item, fig)

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