ids app: implementation of report (WIP)

setup.py

@@ -52,6 +52,7 @@ packages = [
     '%s.plot' % packname,
     '%s.si' % packname,
     '%s.si.ids' % packname,
+    '%s.si.ids.report' % packname,
     '%s.si.ids.targets' % packname]
 
 

src/apps/ids.py

@@ -2,6 +2,7 @@
 # GPLv3
 #
 # The Developers, 21st Century
+import os.path as op
 import sys
 import logging
 from argparse import ArgumentParser
@@ -18,11 +19,15 @@ def d2u(d):
 
 
 subcommand_descriptions = {
-    'go':          'run ids source inversion'
+    'init':        'create ids folder structure and example config',
+    'go':          'run ids source inversion',
+    'report':      'create inversion result report'
 }
 
 subcommand_usages = {
-    'go':          'go <config_file> [options]'
+    'init':        'init <run_dir> [options]',
+    'go':          'go <config_file> [options]',
+    'report':      'report <run_dir> [options]'
 }
 
 subcommands = subcommand_descriptions.keys()
@@ -33,11 +38,13 @@ usage = program_name + ''' <subcommand> <arguments> ... [options]
 
 Subcommands:
 
+    init          %(init)s
     go            %(go)s
+    report        %(report)s
 
 To get further help and a list of available options for any subcommand run:
 
-    mole <subcommand> --help
+    ids <subcommand> --help
 
 ''' % d2u(subcommand_descriptions)
 
@@ -94,6 +101,24 @@ def die(message, err=''):
     sys.exit('%s: error: %s \n %s' % (program_name, message, err))
 
 
+def command_init(args):
+    def setup(parser):
+        parser.add_argument(
+            'rundir', help='new run directory to be created/replaced')
+
+        parser.add_argument(
+            '--force', dest='force', action='store_true',
+            help='overwrite existing run directory (if needed). '
+                 'default is %(default)s')
+
+    parser, args = cl_parse('init', args, setup=setup)
+
+    if args.rundir is None:
+        raise IDSError('running directoy needs to be given')
+
+    IDSRunner.init(rundir=args.rundir, force=args.force)
+
+
 def command_go(args):
     def setup(parser):
         parser.add_argument(
@@ -122,6 +147,22 @@ def command_go(args):
     runner.post_process(force=args.force)
 
 
+def command_report(args):
+    def setup(parser):
+        parser.add_argument(
+            'rundir', help='running directory')
+
+    parser, args = cl_parse('report', args, setup=setup)
+
+    rundir = op.abspath(args.rundir)
+
+    logger.info('Load data from {} ...'.format(rundir))
+    runner = IDSRunner.from_rundir(rundir)
+
+    logger.info('Create IDS report (beta) ...')
+    runner.report()
+
+
 def main(args=None):
     if args is None:
         args = sys.argv[1:]

src/plot/ids.py

@@ -14,7 +14,8 @@ from scipy.interpolate import RegularGridInterpolator as scrgi
 from matplotlib import cm, pyplot as plt, patheffects
 from matplotlib.ticker import FuncFormatter
 
-from pyrocko import gmtpy, orthodrome as pod, moment_tensor as pmt, pile
+from pyrocko import gmtpy, orthodrome as pod, moment_tensor as pmt
+from pyrocko.pile import make_pile
 from pyrocko.guts import Object
 from pyrocko.plot import (
     mpl_color, mpl_graph_color, mpl_init, mpl_margins, mpl_papersize,
@@ -1780,7 +1781,13 @@ def rupture_movie(
 
 
 def trace_view(
-        source, waveform_paths='', scale_mode='common', filename=None, dpi=None):
+        source,
+        pile=None,
+        waveform_paths='',
+        scale_mode='common',
+        filename=None,
+        dpi=None):
+
     mpl_init(fontsize=12)
 
     try:
@@ -1795,12 +1802,18 @@ def trace_view(
         raise ValueError('scale_mode must be chosen out of {}'.format(
             scale_modes))
 
-    p = pile.make_pile(paths=waveform_paths)
+    if pile is not None and waveform_paths:
+        raise ValueError('argmuent collision: "p" and "waveform_paths" are '
+                         'mutually exclusive')
+
+    if waveform_paths:
+        pile = make_pile(paths=waveform_paths)
 
-    channels = num.unique([tr.channel[-1] for tr in p.iter_traces()])
+    channels = num.unique([tr.channel[-1] for tr in pile.iter_traces()])
     channels = num.sort(channels)
 
-    stations = num.unique(['.'.join(tr.nslc_id[:2]) for tr in p.iter_traces()])
+    stations = num.unique([
+        '.'.join(tr.nslc_id[:2]) for tr in pile.iter_traces()])
 
     n_channels = len(channels)
     n_stations = len(stations)
@@ -1821,7 +1834,7 @@ def trace_view(
     scaler = AutoScaler(mode='symmetric', approx_ticks=3.0, space=0.1)
 
     for i, traces in enumerate(
-            p.chopper_grouped(gather=lambda tr: tr.nslc_id[:2])):
+            pile.chopper_grouped(gather=lambda tr: tr.nslc_id[:2])):
         if not traces:
             continue
 

src/scaling.py

@@ -0,0 +1,76 @@
+# GPLv3
+#
+# The Developers, 21st Century
+import logging
+
+import numpy as num
+
+logger = logging.getLogger('ewrica.scaling')
+
+km = 1e3
+
+
+def length_blaser(magnitude, rake):
+    '''Scaling relation for rupture length following  Blaser et al., 2010
+
+    The standard deviation is taken into account leading to slightly different
+    results.
+
+    :param magnitude: Moment magnitude Mw
+    :type magnitude: float
+    :param rake: rake angle in degrees on rupture plane
+    :type rake: float
+
+    :returns: rupture length in m
+    :rtypev float
+    '''
+    rake = (rake % 360.) - 180.
+
+    if rake < 135 and rake > 45:
+        ls = 0.18
+        lcov = num.matrix([[26.14, -3.67], [-3.67, 0.52]]) * 1e-5
+        la, lb = num.random.multivariate_normal([-2.37, 0.57], lcov)
+
+    elif rake > -135 and rake < -45:
+        ls = 0.18
+        lcov = num.matrix([[222.24, -32.34], [-32.34, 4.75]]) * 1e-5
+        la, lb = num.random.multivariate_normal([-1.91, 0.52], lcov)
+
+    else:
+        ls = 0.18
+        lcov = num.matrix([[12.37, -1.94], [-1.94, 0.31]]) * 1e-5
+        la, lb = num.random.multivariate_normal([-2.69, 0.64], lcov)
+
+    return 10**(num.random.normal(magnitude * lb + la, ls**2)) * km
+
+
+def width_blaser(magnitude, rake):
+    '''Scaling relation for rupture width following  Blaser et al., 2010
+
+    The standard deviation is taken into account leading to slightly different
+    results.
+
+    :param magnitude: Moment magnitude Mw
+    :type magnitude: float
+    :param rake: rake angle in degrees on rupture plane
+    :type rake: float
+
+    :returns: rupture width in m
+    :rtypev float
+    '''
+    if rake < 135 and rake > 45:
+        ws = 0.17
+        wcov = num.matrix([[27.47, -3.77], [-3.77, 0.52]]) * 1e-5
+        wa, wb = num.random.multivariate_normal([-1.86, 0.46], wcov)
+
+    elif rake > -135 and rake < -45:
+        ws = 0.16
+        wcov = num.matrix([[264.18, -42.02], [-42.02, 6.73]]) * 1e-5
+        wa, wb = num.random.multivariate_normal([-1.20, 0.36], wcov)
+
+    else:
+        ws = 0.15
+        wcov = num.matrix([[13.48, -2.18], [-2.18, 0.36]]) * 1e-5
+        wa, wb = num.random.multivariate_normal([-1.12, 0.33], wcov)
+
+    return 10**(num.random.normal(magnitude * wb + wa, ws**2)) * km

src/si/ids/core.py

@@ -11,11 +11,14 @@ 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
-from ewrica.io.ids import load_ids_source, load_obssyn
+from ewrica.sources import dump_sources, EwricaIDSSource, load_one_source
+from ewrica.io.ids import (
+    load_ids_source, load_obssyn, load_obssyn_gnss)
 from .config import IDSConfigFull
 from .dataset import Dataset
 from .targets import WaveformTargetGroup, GNSSTargetGroup, InSARTargetGroup
@@ -39,6 +42,8 @@ class IDSRunner(Object):
     config = IDSConfigFull.T(
         default=IDSConfigFull())
 
+    source = EwricaIDSSource.T(optional=True)
+
     keep_tmp = Bool.T(default=False)
 
     def __init__(self, *args, **kwargs):
@@ -94,10 +99,44 @@ class IDSRunner(Object):
     def tracefile(self):
         return op.join(self.resultdir, 'traces.mseed')
 
+    @property
+    def gnssobs_file(self):
+        return op.join(self.resultdir, 'gnss_obs.yaml')
+
+    @property
+    def gnsssyn_file(self):
+        return op.join(self.resultdir, 'gnss_syn.yaml')
+
+    @property
+    def insarfile(self):
+        return op.join(self.resultdir, 'insar.yaml')
+
     @property
     def sourcefile(self):
         return op.join(self.resultdir, 'idssource.yaml')
 
+    @classmethod
+    def from_rundir(cls, rundir):
+        if not op.exists(rundir):
+            raise IDSError('Running directory {} not found'.format(rundir))
+
+        runner = cls()
+        runner._tempdir = rundir
+
+        if not op.exists(op.join(runner.inputdir)):
+            raise IDSError(
+                'Input directory {} not found'.format(runner.inputdir))
+
+        if not op.exists(op.join(runner.rundir)):
+            raise IDSError(
+                'Running directory {} not found'.format(runner.rundir))
+
+        if not op.exists(op.join(runner.resultdir)):
+            raise IDSError(
+                'Result directory {} not found'.format(runner.resultdir))
+
+        return runner
+
     def get_dataset(self):
         if self._dataset is None:
             wf_targets = []
@@ -122,6 +161,22 @@ class IDSRunner(Object):
 
         return self._dataset
 
+    @classmethod
+    def init(cls, rundir, force=False):
+        if op.exists(rundir) and not force:
+            raise FileExistsError(
+                'Given rundir {} already exists. Consider the '
+                '"force" option'.format(rundir))
+        elif op.exists(rundir):
+            shutil.rmtree(rundir)
+
+        os.makedirs(rundir)
+        os.makedirs(op.join(rundir, 'config'))
+
+        runner = cls()
+        runner.config.dump(
+            filename=op.join(rundir, 'config', 'idsconfig_example.conf'))
+
     def _prepare_ids_dir(self):
         dirs = \
             [self.waveformdir, self.gnssdir, self.insardir, self.faultdir] + \
@@ -260,19 +315,29 @@ in the directory %s'''.lstrip() % (
             result_dir = self.resultdir
             os.makedirs(result_dir)
 
+            # Dump source
             dump_sources(
                 [source],
                 filename=self.sourcefile)
 
-            traces = load_obssyn(
-                config_fn=op.join(self.tempdir, self.config.inputfn))
-
+            # 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)
 
+            # Dump observed/modelled GNSS data (if available)
+            try:
+                gnss_obs, gnss_syn = load_obssyn_gnss(config_fn=config_fn)
+                gnss_obs.dump(filename=self.gnssobs_file)
+                gnss_syn.dump(filename=self.gnsssyn_file)
+
+            except FileNotFoundError:
+                pass
+
             if op.exists(outpath_tmpl) and not force:
-                raise IDSError('Configured outpath is already used. Consider '
-                               'using the "force" keyword argument')
+                raise IDSError('Configured outpath exists already. Consider '
+                               'using the "force" argument')
 
             if op.exists(outpath_tmpl):
                 shutil.rmtree(outpath_tmpl)
@@ -286,4 +351,17 @@ in the directory %s'''.lstrip() % (
 
         return source
 
-    # def report(self):
+    def report(self):
+        from ewrica.si.ids.report import create_report
+
+        source = load_one_source(filename=self.sourcefile)
+        pile = make_pile(paths=[self.tracefile])
+
+        gnss_obs, gnss_syn = None, None
+        if op.exists(self.gnssobs_file):
+            gnss_obs = GNSSCampaign.load(filename=self.gnssobs_file)
+
+        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)

src/si/ids/report/__init__.py

@@ -0,0 +1,4 @@
+# GPLv3
+#
+# The Developers, 21st Century
+from .base import *  # noqa

src/si/ids/report/plot.py

@@ -0,0 +1,108 @@
+# 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 ewrica.sources import load_sources, EwricaIDSSource
+from ewrica.io import ids as idsio
+from ewrica.plot import ids as idsplot
+
+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 Plotter(Object):
+    source = EwricaIDSSource.T()
+
+
+class MapPlotter(Plotter):
+    show_grid = Bool.T(default=True)
+    show_topo = Bool.T(default=False)
+
+    @property
+    def basemap(self):
+        source = self.source
+
+        return idsplot.RuptureMap(
+            source=source,
+            lat=source.lat, lon=source.lon,
+            radius=source.length/2.,
+            show_topo=self.show_topo,
+            show_grid=self.show_grid)
+
+
+class FinalSlipMap(MapPlotter):
+    def plot(self):
+        m = self.basemap
+        m.draw_dislocation()
+        m.draw_nucleation_point()
+        m.draw_top_edge()
+        m.save(op.join(run, 'results', 'final_slip_map.png'))
+
+
+class FinalStaticGNSSMap(MapPlotter):
+    def plot(self):
+        gnss_obs, gnss_mod = idsio.load_obssyn_gnss(
+            config_fn=op.join(run, 'input', 'ids_input.inp'))
+
+        m = self.basemap
+        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'))
+
+
+class FinalSlipView(Plotter):
+    def plot(self):
+        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'))
+
+
+class SourceTimeFunction(Plotter):
+    def plot(self):
+        v = idsplot.RuptureView(source=source)
+        v.draw_source_dynamics('stf')
+        v.save(op.join(run, 'results', 'stf.png'))
+
+
+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')
+
+
+plot_classes = [
+    FinalSlipMap,
+    FinalStaticGNSSMap,
+    FinalSlipView,
+    SourceTimeFunction,
+    TracePlot]
+
+
+def create_report_plots():
+    pass
+