Browse Source

ids: WIP

- report command WIP
	- Plotting WIP
	- harvest command WIP
mole_dev
mmetz 3 months ago
parent
commit
f995d74aad
  1. 33
      src/apps/ids.py
  2. 172
      src/io/ids.py
  3. 19
      src/plot/ids.py
  4. 4
      src/si/ids/config.py
  5. 69
      src/si/ids/core.py
  6. 285
      src/si/ids/report/plot.py

33
src/apps/ids.py

@ -2,6 +2,7 @@
# GPLv3
#
# The Developers, 21st Century
import os
import os.path as op
import sys
import logging
@ -21,12 +22,14 @@ def d2u(d):
subcommand_descriptions = {
'init': 'create ids folder structure and example config',
'go': 'run ids source inversion',
'harvest': 'harvest all inversion results (needed to create report)',
'report': 'create inversion result report'
}
subcommand_usages = {
'init': 'init <run_dir> [options]',
'go': 'go <config_file> [options]',
'harvest': 'harvest <run_dir> [options]',
'report': 'report <run_dir> [options]'
}
@ -40,6 +43,7 @@ Subcommands:
init %(init)s
go %(go)s
harvest %(harvest)s
report %(report)s
To get further help and a list of available options for any subcommand run:
@ -143,8 +147,29 @@ def command_go(args):
logger.info('Call IDS backend ...')
runner.run()
logger.info('IDS calculation done. Post process results ...')
runner.post_process(force=args.force)
logger.info('IDS calculation done. Harvest results ...')
runner.harvest(force=args.force)
def command_harvest(args):
def setup(parser):
parser.add_argument(
'rundir', help='running directory')
parser.add_argument(
'--force', dest='force', action='store_true',
help='overwrite existing run directories (if needed). '
'default is %(default)s')
parser, args = cl_parse('harvest', args, setup=setup)
rundir = op.abspath(args.rundir)
logger.info('Load data from {} ...'.format(rundir))
runner = IDSRunner.from_rundir(rundir)
logger.info('Harvest results ...')
runner.harvest(force=args.force)
def command_report(args):
@ -159,8 +184,10 @@ def command_report(args):
logger.info('Load data from {} ...'.format(rundir))
runner = IDSRunner.from_rundir(rundir)
cwd = os.getcwd()
logger.info('Create IDS report (beta) ...')
runner.report()
runner.report(cwd=cwd)
def main(args=None):

172
src/io/ids.py

@ -9,9 +9,10 @@ import numpy as num
from pyrocko import orthodrome as pod, util as putil, moment_tensor as pmt
from pyrocko.cake import earthradius
from pyrocko.guts import Float, Int, List, Object, StringChoice
from pyrocko.guts import Float, Int, List, Object, StringChoice, Tuple
from pyrocko.guts_array import Array
from pyrocko.model import GNSSCampaign, GNSSComponent, GNSSStation
from pyrocko.trace import Trace
from pyrocko.gf.seismosizer import map_anchor
from ewrica.sources import (EwricaIDSSource, EwricaIDSPatch, EwricaIDSSTF,
@ -27,6 +28,68 @@ km = 1e3
d2r = num.pi / 180.
class WaveformResult(Object):
nslc = Tuple.T(optional=True)
observed = Trace.T(optional=True)
synthetic = Trace.T(optional=True)
misfit = Float.T(optional=True)
class WaveformResultList(Object):
results = List.T(
WaveformResult.T(),
default=[],
optional=True)
@classmethod
def primitive_merge(cls, lists):
result = cls()
for lst in lists:
result.results += lst.results
return result
def iter_traces(self):
for res in self.results:
for tr in (res.observed, res.synthetic):
yield tr
def grouped_results(self, gather):
sort = list(map(gather, self.results))
if all(map(lambda x: isinstance(x, bool), sort)):
idx = num.arange(0, len(sort))[sort]
res = [self.results[i] for i in idx]
yield WaveformResultList(results=res)
else:
sort_set = set(sort)
for it in sort_set:
idx = [i for i, s in enumerate(sort) if s == it]
res = [self.results[i] for i in idx]
yield WaveformResultList(results=res)
def grouped_traces(self, gather):
traces = list(self.iter_traces())
sort = list(map(gather, traces))
if all(map(lambda x: isinstance(x, bool), sort)):
idx = num.arange(0, len(sort))[sort]
yield [traces[i] for i in idx]
else:
sort_set = set(sort)
for it in sort_set:
idx = [i for i, s in enumerate(sort) if s == it]
yield [traces[i] for i in idx]
snapshot_parameters = [
'Latitude[deg]',
'Longitude[deg]',
@ -101,6 +164,15 @@ obssyn_gnss_parameters = [
'Up_Syn']
wfmisfitvar_parameters = [
'Station',
'Lat[deg]',
'Lon[deg]',
'nrv_vel_e',
'nrv_vel_n',
'nrv_vel_z']
class IDSUnpackError(Exception):
pass
@ -143,6 +215,10 @@ def _validate_obssyn_gnss_file(fn):
_validate_ids_file(fn, obssyn_gnss_parameters)
def _validate_wf_misfit_file(fn):
_validate_ids_file(fn, wfmisfitvar_parameters)
def load_snapshot_file(fn):
'''Load all information from IDSSnapshot or SnapshotR file
@ -779,17 +855,46 @@ def load_ids_source(config_fn):
return source
def load_wfmisfit_file(fn):
'''Load waveform variance misfit file
:param fn: filename of the WFMisfitvar.dat file
:type fn: str
:returns: station nsl codes, lats, lons, misfit variances of the velocity
in E, N and Z
:rtype: str, float, float, float, float, float
'''
_validate_wf_misfit_file(fn)
stations = num.genfromtxt(
fn,
unpack=True,
skip_header=1,
usecols=(0),
dtype=str)
lat, lon, mf_e, mf_n, mf_z = num.genfromtxt(
fn,
unpack=True,
skip_header=1,
usecols=(1, 2, 3, 4, 5))
return stations, lat, lon, mf_e, mf_n, mf_z
def load_obssyn_file(fn, config_fn):
'''Load one ObsSyn.dat file
:param fn: filename of the ids config file
:param fn: filename of the obssyn file
:type fn: str
:param config_fn: filename of the ids config file
:type config_fn: str
:returns: pyrocko trace
:rtype: :py:class:`pyrocko.trace.Trace`
'''
from pyrocko.trace import Trace
try:
_validate_obssyn_nearfield_file(fn)
nearfield = True
@ -822,28 +927,32 @@ def load_obssyn_file(fn, config_fn):
deltat = num.mean(num.diff(time_h))
def generate_traces(data, channels, obs_syn, tmin):
return [
Trace(
def generate_traces(obs, syn, channel, tmin):
return WaveformResult(
observed=Trace(
tmin=tmin,
deltat=deltat,
ydata=d,
ydata=obs,
network=network,
station=station,
location=obs_syn,
channel=c)
for d, c in zip(data, channels)]
location=location,
channel=channel),
synthetic=Trace(
tmin=tmin,
deltat=deltat,
ydata=syn,
network=network,
station=station,
location=location,
channel=channel),
nslc=(network, station, location, channel))
traces = generate_traces(
data=[obs_e, obs_n], channels=['E', 'N'], obs_syn='OB', tmin=tmin_h)
traces += generate_traces(
data=[obs_z], channels=['Z'], obs_syn='OB', tmin=tmin_z)
traces += generate_traces(
data=[syn_e, syn_n], channels=['E', 'N'], obs_syn='SY', tmin=tmin_h)
traces += generate_traces(
data=[syn_z], channels=['Z'], obs_syn='SY', tmin=tmin_z)
results = []
results.append(generate_traces(obs_e, syn_e, 'E', tmin_h))
results.append(generate_traces(obs_n, syn_n, 'N', tmin_h))
results.append(generate_traces(obs_z, syn_z, 'Z', tmin_z))
return traces
return WaveformResultList(results=results)
def load_obssyn(config_fn):
@ -858,13 +967,28 @@ def load_obssyn(config_fn):
config = load_ids_config_file(config_fn)
rundir = op.join(op.dirname(config_fn), config['rundir'])
wfmisfit_fn = op.join(rundir, config['wfvariance_file'])
stations, _, _, mf_e, mf_n, mf_z = load_wfmisfit_file(wfmisfit_fn)
misfits = dict(
E=mf_e,
N=mf_n,
Z=mf_z)
traces = []
results = []
for fn in os.listdir(op.join(rundir, config['o2sdir'])):
traces += load_obssyn_file(
op.join(rundir, config['o2sdir'], fn), config_fn)
results.append(load_obssyn_file(
op.join(rundir, config['o2sdir'], fn), config_fn))
waveform_result = WaveformResultList.primitive_merge(results)
for result in waveform_result.results:
idx = num.where(stations == '.'.join(result.nslc[:3]))[0]
if num.any(idx):
result.misfit = misfits[result.nslc[3]][num.asscalar(idx)]
return traces
return waveform_result
def load_obssyn_gnss(config_fn):

19
src/plot/ids.py

@ -1786,7 +1786,8 @@ def trace_view(
waveform_paths='',
scale_mode='common',
filename=None,
dpi=None):
dpi=None,
show_fig=False):
mpl_init(fontsize=12)
@ -1803,7 +1804,7 @@ def trace_view(
scale_modes))
if pile is not None and waveform_paths:
raise ValueError('argmuent collision: "p" and "waveform_paths" are '
raise ValueError('argument collision: "pile" and "waveform_paths" are '
'mutually exclusive')
if waveform_paths:
@ -1818,8 +1819,8 @@ def trace_view(
n_channels = len(channels)
n_stations = len(stations)
fig = plt.figure(figsize=mpl_papersize(paper='a4', orientation='portrait'))
fig.set_figheight(50. * (n_stations + 2) / 72)
figsize = mpl_papersize(paper='a4', orientation='portrait')
fig = plt.figure(figsize=(figsize[0], 50. * (n_stations + 2) / 72))
mpl_margins(
fig,
@ -1847,6 +1848,11 @@ def trace_view(
tr_syn = tr
elif tr.location == 'OB':
tr_obs = tr
elif tr.meta is not None:
if tr.meta['obs_syn'] == 'SYN':
tr_syn = tr
elif tr.meta['obs_syn'] == 'OBS':
tr_obs = tr
sharex = None
sharey = None
@ -1940,9 +1946,12 @@ def trace_view(
plt.clf()
plt.close()
else:
elif show_fig:
plt.show()
else:
return fig
__all__ = [
'make_colormap',

4
src/si/ids/config.py

@ -61,6 +61,10 @@ class RunConfig(Object):
def outpath_template(self, outpath_template):
self.outpath_template__ = outpath_template
@property
def run(self):
return op.basename(self.outpath_template)
class EngineConfig(Object):
ids_store_superdirs = List.T(String.T(), default=['.'])

69
src/si/ids/core.py

@ -11,14 +11,12 @@ from subprocess import Popen, PIPE
import os
import os.path as op
from pyrocko.pile import make_pile
from pyrocko.guts import Object, Bool
from pyrocko.io import save
from pyrocko.model.gnss import GNSSCampaign
from ewrica.sources import dump_sources, EwricaIDSSource, load_one_source
from ewrica.io.ids import (
load_ids_source, load_obssyn, load_obssyn_gnss)
load_ids_source, load_obssyn, load_obssyn_gnss, WaveformResultList)
from .config import IDSConfigFull
from .dataset import Dataset
from .targets import WaveformTargetGroup, GNSSTargetGroup, InSARTargetGroup
@ -95,9 +93,14 @@ class IDSRunner(Object):
def resultdir(self):
return op.join(self.tempdir, 'results')
@property
def reportdir(self):
return 'report'
@property
def tracefile(self):
return op.join(self.resultdir, 'traces.mseed')
return op.join(self.resultdir, 'waveform_fits.yaml')
# return op.join(self.resultdir, 'traces.mseed')
@property
def gnssobs_file(self):
@ -115,6 +118,14 @@ class IDSRunner(Object):
def sourcefile(self):
return op.join(self.resultdir, 'idssource.yaml')
@property
def configfile(self):
return op.join(self.tempdir, self.config.inputdir, 'ids_config.yaml')
@property
def _input_fn(self):
return op.join(self.tempdir, self.config.inputfn)
@classmethod
def from_rundir(cls, rundir):
if not op.exists(rundir):
@ -122,6 +133,8 @@ class IDSRunner(Object):
runner = cls()
runner._tempdir = rundir
runner.config = IDSConfigFull()
runner.config = IDSConfigFull.load(filename=runner.configfile)
if not op.exists(op.join(runner.inputdir)):
raise IDSError(
@ -200,8 +213,9 @@ class IDSRunner(Object):
def run(self):
config = self.config
input_fn = op.join(self.tempdir, config.inputfn)
self._input_fn = input_fn
# input_fn = op.join(self.tempdir, config.inputfn)
# self._input_fn = input_fn
input_fn = self._input_fn
with open(input_fn, 'wb') as f:
input_str = config.string_for_config(self.get_dataset())
@ -305,14 +319,28 @@ in the directory %s'''.lstrip() % (
os.chdir(old_wd)
def post_process(self, force=False):
def harvest(self, force=False):
source = load_ids_source(self._input_fn)
source.name = self.config.dataset_config.event_name
outpath_tmpl = self.config.run_config.outpath_template
if outpath_tmpl is not None:
logger.info('Move run dir to {} ...'.format(outpath_tmpl))
if op.exists(outpath_tmpl) and not force:
raise IDSError('Configured outpath exists already. Consider '
'using the "force" argument')
if self.tempdir == outpath_tmpl:
self._tempdir = ''
shutil.copytree(outpath_tmpl, self.tempdir, dirs_exist_ok=True)
result_dir = self.resultdir
if op.exists(result_dir):
shutil.rmtree(result_dir)
os.makedirs(result_dir)
# Dump source
@ -322,9 +350,11 @@ in the directory %s'''.lstrip() % (
# Dump observed/modelled traces
config_fn = op.join(self.tempdir, self.config.inputfn)
traces = load_obssyn(config_fn=config_fn)
save(
traces, filename_template=self.tracefile)
waveform_results = load_obssyn(config_fn=config_fn)
waveform_results.dump(filename=self.tracefile)
# save(
# traces, filename_template=self.tracefile)
# Dump observed/modelled GNSS data (if available)
try:
@ -335,10 +365,6 @@ in the directory %s'''.lstrip() % (
except FileNotFoundError:
pass
if op.exists(outpath_tmpl) and not force:
raise IDSError('Configured outpath exists already. Consider '
'using the "force" argument')
if op.exists(outpath_tmpl):
shutil.rmtree(outpath_tmpl)
@ -351,11 +377,11 @@ in the directory %s'''.lstrip() % (
return source
def report(self):
def report(self, cwd='.'):
from ewrica.si.ids.report import create_report
source = load_one_source(filename=self.sourcefile)
pile = make_pile(paths=[self.tracefile])
waveform_result = WaveformResultList.load(filename=self.tracefile)
gnss_obs, gnss_syn = None, None
if op.exists(self.gnssobs_file):
@ -364,4 +390,13 @@ in the directory %s'''.lstrip() % (
if op.exists(self.gnsssyn_file):
gnss_syn = GNSSCampaign.load(filename=self.gnsssyn_file)
create_report(source, pile, gnss_obs=gnss_obs, gnss_syn=gnss_syn)
reportdir = op.join(cwd, self.reportdir)
run = self.config.run_config.run
create_report(
reportdir,
run,
source,
waveform_result,
gnss_obs=gnss_obs,
gnss_syn=gnss_syn)

285
src/si/ids/report/plot.py

@ -1,19 +1,13 @@
# GPLv3
#
# The Developers, 21st Century
import os
import os.path as op
import numpy as num
import matplotlib.pyplot as plt
from pyrocko import model
from pyrocko.guts import Bool, Object
from pyrocko.plot.dynamic_rupture import make_colormap
from pyrocko.plot import mpl_init, mpl_papersize
from pyrocko.guts import Bool, Object, String
from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize, mpl_graph_color
from pyrocko.util import tts
from ewrica.sources import load_sources, EwricaIDSSource
from ewrica.io import ids as idsio
from ewrica.plot import ids as idsplot
mpl_init(fontsize=15)
@ -21,88 +15,277 @@ mpl_init(fontsize=15)
km = 1e3
runs = [d for d in os.listdir('.') if op.isdir(d) and 'curved' not in d]
for run in runs:
source = load_sources(filename=op.join(run, 'results', 'idssource.yaml'))[0]
print(run, source.misfit)
class PlotGroup(Object):
title = String.T(default='')
description = String.T(default='')
class PlotItem(Object):
title = String.T(default='')
description = String.T(default='')
name = String.T(default='')
class Plot(Object):
name = String.T(default='')
def make(self, source, *args, **kwargs):
return self, PlotGroup(), self.draw_figures(source, *args, **kwargs)
def draw_figures(self, source, *args, **kwargs):
raise NotImplementedError('Method not implemented')
class Plotter(Object):
source = EwricaIDSSource.T()
class MPLPlot(Plot):
pass
class ViewPlot(Plot):
pass
class MapPlotter(Plotter):
class GMTPlot(Plot):
show_grid = Bool.T(default=True)
show_topo = Bool.T(default=False)
@property
def basemap(self):
source = self.source
def basemap(self, source):
return idsplot.RuptureMap(
source=source,
lat=source.lat, lon=source.lon,
radius=source.length/2.,
radius=source.length,
show_topo=self.show_topo,
show_grid=self.show_grid)
class FinalSlipMap(MapPlotter):
def plot(self):
m = self.basemap
class FinalSlipMap(GMTPlot):
name = 'slipmap'
def make(self, source, *args, **kwargs):
return (
self,
PlotGroup(
title='Final slip map',
description=''),
self.draw_figures(source, *args, **kwargs))
def draw_figures(self, source, *args, **kwargs):
m = self.basemap(source)
m.draw_dislocation()
m.draw_nucleation_point()
m.draw_top_edge()
m.save(op.join(run, 'results', 'final_slip_map.png'))
item = PlotItem(
title='best slip distribution map',
description='',
name='slip_map')
class FinalStaticGNSSMap(MapPlotter):
def plot(self):
gnss_obs, gnss_mod = idsio.load_obssyn_gnss(
config_fn=op.join(run, 'input', 'ids_input.inp'))
yield item, m
m = self.basemap
class FinalStaticGNSSMap(GMTPlot):
name = 'staticgnss_map'
def make(self, source, *args, **kwargs):
return (
self,
PlotGroup(
title='Static GNSS Fit',
description=''),
self.draw_figures(source, *args, **kwargs))
def draw_figures(self, source, gnss_obs, gnss_syn):
m = self.basemap(source)
m.draw_dislocation()
m.draw_nucleation_point()
m.draw_top_edge()
m.draw_gnss_data_fit(
campaign_obs=gnss_obs, campaign_syn=gnss_mod, vertical=False)
m.save(op.join(run, 'results', 'gnss_fit_horizontal_map.png'))
campaign_obs=gnss_obs, campaign_syn=gnss_syn, vertical=False)
item = PlotItem(
title='static gnss fit map',
description='',
name='gnss_fit')
yield item, m
class FinalSlipView(ViewPlot):
name = 'slipview'
def make(self, source, *args, **kwargs):
return (
self,
PlotGroup(
title='Final Slip View',
description=''),
self.draw_figures(source, *args, **kwargs))
class FinalSlipView(Plotter):
def plot(self):
def draw_figures(self, source, *args, **kwargs):
v = idsplot.RuptureView(source=source)
v.draw_dislocation()
v.draw_nucleation_point()
v.draw_dislocation_contour()
v.save(op.join(run, 'results', 'final_slip_view.png'))
item = PlotItem(
title='best slip distribution view',
description='',
name='slip_view')
class SourceTimeFunction(Plotter):
def plot(self):
yield item, v
class SourceTimeFunction(ViewPlot):
name = 'stf'
def make(self, source, *args, **kwargs):
return (
self,
PlotGroup(
title='Source Time Function',
description=''),
self.draw_figures(source, *args, **kwargs))
def draw_figures(self, source, *args, **kwargs):
v = idsplot.RuptureView(source=source)
v.draw_source_dynamics('stf')
v.save(op.join(run, 'results', 'stf.png'))
item = PlotItem(
title='source time function of best model',
description='',
name='stf')
class TracePlot(Plotter):
def plot(self):
idsplot.trace_view(
source=source,
waveform_paths=[op.join(run, 'results', 'traces.mseed')],
filename=op.join(run, 'results', 'traces.png'),
scale_mode='individual')
yield item, v
class TracePlot(MPLPlot):
name = 'waveforms'
def make(self, source, *args, **kwargs):
return (
self,
PlotGroup(
title='Waveform Fit',
description='Black - observed, red - modelled'),
self.draw_figures(source, *args, **kwargs))
def draw_figures(self, source, waveform_result, *args, **kwargs):
# TODO write single component single station plot
try:
plt.rcParams['ytick.right'] = True
plt.rcParams['ytick.labelright'] = True
plt.rcParams['ytick.left'] = False
plt.rcParams['ytick.labelleft'] = False
except KeyError:
plt.rcParams['ytick.right'] = True
plt.rcParams['ytick.left'] = False
channels = num.unique([
tr.channel[-1] for tr in waveform_result.iter_traces()])
channels = num.sort(channels)
axes = []
for i, res_sta in enumerate(waveform_result.grouped_results(
gather=lambda res: res.nslc[:3])):
if not res_sta:
continue
for ic, res_cha in enumerate(res_sta.grouped_results(
gather=lambda res: res.nslc)):
res_cha = res_cha.results[0]
c = res_cha.nslc[3]
tr_obs = res_cha.observed
tr_syn = res_cha.synthetic
fig = plt.figure(
figsize=mpl_papersize(paper='a6', orientation='landscape'))
mpl_margins(
fig,
top=4., bottom=40.,
left=4., right=75.,
units='point')
ax = fig.add_subplot(111)
line_obs, = ax.plot(
tr_obs.get_xdata() - source.time, tr_obs.get_ydata(),
c='dimgrey')
line_syn, = ax.plot(
tr_syn.get_xdata() - source.time, tr_syn.get_ydata(),
c=mpl_graph_color(0))
ax.tick_params(direction="in", width=1.5, length=6)
if any(tr_obs.get_xdata() < source.time):
ax.axvline(
x=0.,
color='darkgrey')
ax.axhline(
y=0.,
color='lightgrey',
zorder=0.)
text_kwargs = dict(
x=0.95, y=0.15,
s='{}\nmisfit = {} m/s'.format(
'.'.join(tr_syn.nslc_id[:2] + (c,)),
res_cha.misfit),
transform=ax.transAxes,
horizontalalignment='right',
verticalalignment='center')
try:
ax.text(
fontfamily='monospace',
**text_kwargs)
except TypeError:
from matplotlib.font_manager import FontProperties
font0 = FontProperties()
font0.set_family('monospace')
ax.text(
fontproperties=font0,
**text_kwargs)
line_obs.set_label('observed')
line_syn.set_label('modelled')
plt.legend(
ncol=2,
bbox_to_anchor=(-0.5, -0.9),
borderaxespad=0.,
loc='center')
ax.set_xlabel('seconds after {}'.format(tts(source.time)))
ax.set_ylabel('velocity [m/s]')
ax.yaxis.set_label_position("right")
axes.append(ax)
item = PlotItem(
title='trace fit observed to modelled ',
description='',
name='trace_fit_{}'.format(
'.'.join(tr_obs.nslc_id[:2] + (c,))))
yield item, fig
plt.cla()
plt.clf()
plt.close()
# TODO Station Distribution
plot_classes = [
FinalSlipMap,
FinalStaticGNSSMap,
FinalSlipView,
SourceTimeFunction,
TracePlot]
TracePlot,
SourceTimeFunction]
def create_report_plots():
pass
plot_classes_gnss = [
FinalStaticGNSSMap]
Loading…
Cancel
Save