Browse Source

ids wrapper: bugfixing and testing vs old version

pull/4/head
mmetz 2 years ago
parent
commit
2117ac48a6
  1. 951
      src/si/ids/config.py
  2. 503
      src/si/ids/core.py
  3. 799
      src/si/ids/dataset.py
  4. 207
      src/si/ids/targets/gnss.py
  5. 518
      src/si/ids/targets/waveform.py
  6. 36
      test/io/test_ids.py

951
src/si/ids/config.py

@ -1,473 +1,478 @@
# GPLv3
#
# The Developers, 21st Century
import logging
import os
import os.path as op
import numpy as num
from pyrocko import util as putil
from pyrocko.guts import Float, List, String, Dict, Int, Object
from .targets import (
TargetGroupConfig, WaveformTargetGroupConfig, GNSSTargetGroupConfig,
InSARTargetGroupConfig)
from .dataset import DatasetConfig
logger = logging.getLogger('ewrica.si.ids.config')
km = 1e3
class RunConfig(Object):
n_iterations = Int.T(default=100)
outpath_template__ = String.T(optional=True)
deltat_snapshots = Float.T(default=5.)
t_max_snapshots = Float.T(default=50.)
_expand_paths = TargetGroupConfig._expand_paths
def __init__(self, *args, parent=None, **kwargs):
Object.__init__(self, *args, **kwargs)
self._parent = parent
def bind_parent(self, parent):
self._parent = parent
def unbind_parent(self):
self._parent = None
@staticmethod
def example():
return RunConfig()
@property
def outpath_template(self):
return self._expand_paths(self.outpath_template__)
@outpath_template.setter
def outpath_template(self, outpath_template):
self.outpath_template__ = outpath_template
class EngineConfig(Object):
ids_store_superdirs = List.T(String.T(), default=['.'])
pyrocko_store_superdirs = List.T(String.T(), default=['.'])
ids_waveform_store_id = String.T(default='.')
ids_geodetic_store_id = String.T(default='.')
ids_geodetic_files = Dict.T(
String.T(), String.T(),
default=dict(Z='uz', R='ur', T='ut'))
ids_synthetic_bandpass_filter = Dict.T(
String.T(), Float.T(),
default=dict(order=3., f_min=0.01, f_max=0.5))
@staticmethod
def example():
return EngineConfig(
ids_store_superdirs=['../ids_store_superdirs'],
pyrocko_store_superdirs=['../pyrocko_store_superdirs'],
ids_waveform_store_id='ids_waveform_store',
ids_geodetic_store_id='ids_geodetic_store',
ids_synthetic_bandpass_filter=dict(
order=3,
f_min=0.01,
f_max=0.5))
@property
def ids_geodetic_store(self):
for d in self.ids_store_superdirs:
store = op.join(d, self.ids_geodetic_store_id)
if op.exists(store):
if not store.endswith(os.sep):
store += os.sep
return store
@property
def ids_waveform_store(self):
for d in self.ids_store_superdirs:
store = op.join(d, self.ids_waveform_store_id)
if op.exists(store):
if not store.endswith(os.sep):
store += os.sep
return store
class IDSConfig(Object):
ids_version = String.T(default='2020')
path_prefix = String.T(
default=op.expanduser('~'),
help='Location of parent directory. All other dataset paths can be '
'given as relative paths from "path_prefix"',
optional=True)
dataset_config = DatasetConfig.T(default=DatasetConfig())
waveform_config = List.T(
WaveformTargetGroupConfig.T(),
default=[WaveformTargetGroupConfig()])
gnss_config = GNSSTargetGroupConfig.T(
default=GNSSTargetGroupConfig(), optional=True)
insar_config = InSARTargetGroupConfig.T(
default=InSARTargetGroupConfig(), optional=True)
run_config = RunConfig.T(default=RunConfig())
engine_config = EngineConfig.T(default=EngineConfig())
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self.bind_parent()
def bind_parent(self):
for conf in (
self.dataset_config,
self.gnss_config,
self.insar_config,
self.run_config):
if conf is None:
continue
conf.bind_parent(self)
for conf in self.waveform_config:
conf.bind_parent(self)
class IDSConfigFull(IDSConfig):
@property
def inputdir(self):
return 'input'
@property
def inputfn(self):
return op.join(self.inputdir, 'ids_input.inp')
@property
def waveformdir(self):
return op.join(self.inputdir, 'waveform')
@property
def gnssdir(self):
return op.join(self.inputdir, 'gnss')
@property
def gnssfile(self):
fn = 'gnss_data.dat'
if self.gnss_config:
fn = self.gnss_config.fn_dat
return op.join(self.gnssdir, fn)
@property
def insardir(self):
return op.join(self.inputdir, 'insar')
@property
def rundir(self):
return 'run'
@staticmethod
def example():
dataset_config = DatasetConfig.example()
waveform_config = [WaveformTargetGroupConfig.example()]
gnss_config = GNSSTargetGroupConfig.example()
insar_config = InSARTargetGroupConfig.example()
engine_config = EngineConfig.example()
run_config = RunConfig.example()
return IDSConfigFull(
path_prefix='..',
dataset_config=dataset_config,
waveform_config=waveform_config,
gnss_config=gnss_config,
insar_config=insar_config,
engine_config=engine_config,
run_config=run_config)
def string_for_config(self, dataset):
event = dataset.event
eng_conf = self.engine_config
np = self.dataset_config.fault_plane[-1]
gnss_target = dataset.gnss_targetgroup
d = dict(
origin_time=putil.tts(event.time, format='%Y %m %d %H %M %S'),
lat=event.lat,
lon=event.lon,
depth=event.depth / km,
mag=event.magnitude,
waveform_dir=op.relpath(
self.waveformdir, start=self.inputdir) + os.sep,
dist_min=dataset.get_waveform_distance_range()[0] / km,
dist_max=dataset.get_waveform_distance_range()[1] / km,
waveform_static_weight=dataset.get_waveform_static_weight(),
gf_dir=eng_conf.ids_waveform_store,
filter_order=eng_conf.ids_synthetic_bandpass_filter['order'],
filter_f_min=eng_conf.ids_synthetic_bandpass_filter['f_min'],
filter_f_max=eng_conf.ids_synthetic_bandpass_filter['f_max'],
gnss_nsta=gnss_target.n_stations if gnss_target else 0,
gnss_dat=op.relpath(
self.gnssfile, start=self.inputdir), # TODO
gnss_weight=dataset.gnss_targetgroup.config.weight, # TODO
insar_ngrid=self.insar_config.n_grids, # TODO
insar_dat=self.insar_config.dat_path, # TODO
insar_weight=self.insar_config.weight, # TODO
geodetic_gf_dir=eng_conf.ids_geodetic_store,
geodetic_gf_z=eng_conf.ids_geodetic_files['Z'],
geodetic_gf_radial=eng_conf.ids_geodetic_files['R'],
geodetic_gf_transvers=eng_conf.ids_geodetic_files['T'],
idisc=1,
subfault='#',
subfault_file="' '",
n_subfaults=0,
ref_subfault=None,
focmec='',
strike=getattr(event.moment_tensor, 'strike{}'.format(np)),
dip=getattr(event.moment_tensor, 'dip{}'.format(np)),
rake=getattr(event.moment_tensor, 'rake{}'.format(np)),
n_iter_max=self.run_config.n_iterations,
run_dir=op.relpath(self.rundir, start=self.inputdir) + os.sep)
t_max = self.run_config.t_max_snapshots
deltat = self.run_config.deltat_snapshots
snaps = ["'Snapshot_{}s.dat' {:<10.1f}{:<10.1f}".format(t, 0.0, t)
for t in num.arange(deltat, t_max + deltat, deltat)]
snaps += [
"'SnapshotR_{}s.dat' {:<10.1f}{:<10.1f}".format(t, t - deltat, t)
for t in num.arange(deltat, t_max + deltat, deltat)]
d['n_snaps'] = len(snaps)
d['snapshot_lines'] = '\n '.join(snaps)
template = '''# Example for the input file of FORTRAN77 program "ids2020" for kinematic earthquake
# source imaging based on the Iterative Deconvolution and Stacking (IDS) Method using
# seismic waveform and static displacement data jointly.
#
# written by Rongjiang Wang
# Helmholtz Centre Potsdam
# GFZ German Research Centre for Geosciences
# e-mail: wang@gfz-potsdam.de
# Phone +49 331 2881209
# Fax +49 331 2881204
#
# Last modified: Potsdam, August, 2020
#
# Reference:
# Zhang and Wang et al. (2014). Automatic imaging of earthquake rupture processes
# by iterative deconvolution and stacking of high-rate GNSS and strong motion
# seismograms, JGR Solid Earth, 199(7), 5633-5650.
#
#################################################################
## ##
## If not specified otherwise, SI Unit System is used overall! ##
## ##
## ##
## Each comment line should start with "#" ##
## and no blank line is allowed! ##
## ##
#################################################################
#
#================================================================================
# EARTHQUAKE LOCATION PARAMETERS
# ------------------------------
# 1. origin time [s] (year, month, day, hour, minute [positive integer number]
# and second [positive or negative real number]), which will be also used as
# the reference time of the output data
# 2. hypocentre (lat[deg], lon[deg], depth[km])
# 3. preliminary estimate of moment magnitude
# Note: the hypocentre location and origin time (at least to minute) should be
# consistent with those given in the "StationInfo.dat" file of the
# waveform data (s. below)
#--------------------------------------------------------------------------------
{origin_time}
{lat} {lon} {depth}
{mag}
#================================================================================
# WAVEFORM DATA (HIGH-RATE-GNSS/STRONG-MOTION SEISMOGRAMS)
# --------------------------------------------------------------------
# 1. folder name of the data
# Note: this folder should include a file 'StationInfo.dat' giving all
# necessary information about the local monitoring network.
# 2. lower and upper thresholds of the epicentral distances [km] (stations
# outside the threshold distances will not be used)
# 3. weight of the static offset data (obtained automatically from the velocity
# waveforms after an empirical baseline correction) relative to the waveform
# data (suggested value 0.0-1.0)
#--------------------------------------------------------------------------------
'{waveform_dir}'
{dist_min} {dist_max}
{waveform_static_weight}
#================================================================================
# GREEN'S FUNCTION DATABASE
# -------------------------
# 1. Green's function database (folder name)
# 2. Parameters of Butterworth bandpass filter applied to the synthetics:
# order, lower and upper cutoff frequencies [Hz]
# computation [Hz]
#--------------------------------------------------------------------------------
'{gf_dir}'
{filter_order} {filter_f_min} {filter_f_max}
#================================================================================
# GNSS STATIC OFFSETS
# -------------------
# 1. number of GNSS sites, file name of the data
# Note: if the number of GNSS sites <= 0, then no GNSS data are available.
# the GNSS data file includes one header line and 8 culomns,
# Lat(deg) Lon(deg) East(m) North(m) Up(m) wf_E wf_N wf_Z
# where wf_E/N/Z are within-dataset relative weighting factor proportional to
# inverse error variance
# 2. weight of the whole GNSS dataset relative to the waveform dataset (suggested
# value 0.0-1.0)
#--------------------------------------------------------------------------------
{gnss_nsta} '{gnss_dat}'
{gnss_weight}
#================================================================================
# InSAR LOS DISPLACEMENTS
# -----------------------
# 1. number of InSAR grids, file name of the data
# Note: if the number of InSAR grids <= 0, then no InSAR data are used.
# the InSAR data file includes one header line and 8 culomns,
# Lat(deg) Lon(deg) LOS(m) wf_los incidence(deg) azimuth(deg)
# where wf_los is the within-dataset relative weighting factor
# proportional to inverse error variance
# 2. weight of the whole InSAR dataset in the misfit variance relative to the
# waveform dataset (suggested value 0.0-1.0)
#--------------------------------------------------------------------------------
{insar_ngrid} '{insar_dat}'
{insar_weight}
#================================================================================
# Geodetic GREEN'S FUNCTION DATABASE
# ----------------------------------
# 1. geodetic Green's function database (folder name)
# 2. 3 file names of the geodetic Green's functions
#--------------------------------------------------------------------------------
'{geodetic_gf_dir}'
'{geodetic_gf_z}' '{geodetic_gf_radial}' '{geodetic_gf_transvers}'
#================================================================================
# FAULT PARAMETERS
# ----------------
# 1. ipatch: switch (0/1) for subfaults input form,
# 0 = read discrete patches from files,
# 1 = automatic discretization of a single rectangular fault
#--------------------------------------------------------------------------------
# idisc
#--------------------------------------------------------------------------------
{idisc}
#--------------------------------------------------------------------------------
# IF(idisc = 0)THEN
#
# subfault paramters
# ==================
# file_name: data file for discrete patches
# number of sub-faults
# iref selection of the patch reference location
#
# Note: the file includes one header line followed by 8 columns of data:
#
# latitude[deg], longitude[deg], depth[km], length[km], width[km],
# strike[deg], dip[deg], rake[deg]
#
# where (lat, lon, dep) is the reference location of the patch depending
# on iref: 1 = upper-left corner, 2 = upper-right corner,
# 3 = lower-left corner, 4 = lower-right corner,
# 5 = central point (s. Fig.)
#
# definitions for a rectangular fault patch
# =========================================
#
# north(x)
# /
# / | strike
# 1-----------------------2
# |\ p . \ W
# :-\ i . \ i
# | \ l . 5 \ d
# :90 \ S . \ t
# |-dip\ . \ h
# : \. ) rake \
# downward(z) 3-----------------------4
# L e n g t h
#--------------------------------------------------------------------------------
# file_name nsf iref
#--------------------------------------------------------------------------------
{subfault} {subfault_file} {n_subfaults} {ref_subfault}
#--------------------------------------------------------------------------------
# ELSE IF(idisc = 1)THEN
#
# 2. focal mechanism: strike, dip and rake angles [deg]
#--------------------------------------------------------------------------------
{focmec} {strike:8.3f} {dip:7.3f} {rake:8.3f}
#================================================================================
# IDS PARAMETERS
# --------------
# 1. maximum numbers of iterations.
# Note: if max. iteration = 0, forward modelling for any dummy data set is
# made using the given rupture model and the filter parameters
# if max. iteration < 0, same as iteration = 0, but without filtering.
#--------------------------------------------------------------------------------
{n_iter_max}
#================================================================================
# OUTPUT
# ------
# 1. folder name for the output files (should exist already)
# 2. sub-folder name for comparisons of observed and modelled waveform data
# (should exist already)
# 3. sub-folder name for results from empirical baseline correction of
# waveform data (this folder should exist already)
# 4. file name of major earthquake parameters
# 5. file name of earthquake source-time-function (moment rate history)
# 6. file name of sub-fault source-time-function (moment rate history)
# 7. file name of rupture propagation
# 8. file name of final slip model
# Note: if forward modelling is selected (s. above), the above two files
# should provide data for the given kinematic rupture model and will
# be overwritten after the IDS inversion
# 9. file name of normalized residual waveform variance
#10. file name of predicted static offsets at seismic waveform data sites
#11. file name of comparison between observed and modelled GNSS data
#12. file name of comparison between observed and modelled InSAR data
#13. number of snapshots of the rupture process
#14. file name of the 1. snapshot for the slip released within the given
# time window (t1[s], t2[s])
#--------------------------------------------------------------------------------
'{run_dir}'
'O2S'
'EBC'
'EQ_MP.dat'
'EQ_STF.dat'
'PT_STF.dat'
'Rup_Front.dat'
'SlipModel.dat'
'WFMisfitvar.dat'
'SMStaticOffset.dat'
'GNSSDatafit.dat'
'InSARDatafit.dat'
{n_snaps}
{snapshot_lines}
#================================end of input====================================
''' # noqa
return template.format(**d).encode('ascii')
__all__ = [
'IDSConfig',
'IDSConfigFull',
'DatasetConfig',
'EngineConfig',
'RunConfig']
# GPLv3
#
# The Developers, 21st Century
import logging
import os
import os.path as op
import numpy as num
from pyrocko import util as putil
from pyrocko.guts import Float, List, String, Dict, Int, Object
from .targets import (
TargetGroupConfig, WaveformTargetGroupConfig, GNSSTargetGroupConfig,
GNSSTargetGroupConfigFull, InSARTargetGroupConfig)
from .dataset import DatasetConfig
logger = logging.getLogger('ewrica.si.ids.config')
km = 1e3
class FilterSetting(Object):
order = Int.T(default=3)
f_min = Float.T(default=0.0)
f_max = Float.T(default=0.5)
class RunConfig(Object):
n_iterations = Int.T(default=100)
outpath_template__ = String.T(optional=True)
deltat_snapshots = Float.T(default=5.)
t_max_snapshots = Float.T(default=50.)
_expand_paths = TargetGroupConfig._expand_paths
def __init__(self, *args, parent=None, **kwargs):
Object.__init__(self, *args, **kwargs)
self._parent = parent
def bind_parent(self, parent):
self._parent = parent
def unbind_parent(self):
self._parent = None
@staticmethod
def example():
return RunConfig()
@property
def outpath_template(self):
return self._expand_paths(self.outpath_template__)
@outpath_template.setter
def outpath_template(self, outpath_template):
self.outpath_template__ = outpath_template
class EngineConfig(Object):
ids_store_superdirs = List.T(String.T(), default=['.'])
pyrocko_store_superdirs = List.T(String.T(), default=['.'])
ids_waveform_store_id = String.T(default='.')
ids_geodetic_store_id = String.T(default='.')
ids_geodetic_files = Dict.T(
String.T(), String.T(),
default=dict(Z='uz', R='ur', T='ut'))
ids_synthetic_bp_filter = FilterSetting.T(default=FilterSetting())
@staticmethod
def example():
return EngineConfig(
ids_store_superdirs=['../ids_store_superdirs'],
pyrocko_store_superdirs=['../pyrocko_store_superdirs'],
ids_waveform_store_id='ids_waveform_store',
ids_geodetic_store_id='ids_geodetic_store')
@property
def ids_geodetic_store(self):
for d in self.ids_store_superdirs:
store = op.join(d, self.ids_geodetic_store_id)
if op.exists(store):
if not store.endswith(os.sep):
store += os.sep
return store
@property
def ids_waveform_store(self):
for d in self.ids_store_superdirs:
store = op.join(d, self.ids_waveform_store_id)
if op.exists(store):
if not store.endswith(os.sep):
store += os.sep
return store
class IDSConfig(Object):
ids_version = String.T(default='2020')
path_prefix = String.T(
default=op.expanduser('~'),
help='Location of parent directory. All other dataset paths can be '
'given as relative paths from "path_prefix"',
optional=True)
dataset_config = DatasetConfig.T(default=DatasetConfig())
waveform_config = List.T(
WaveformTargetGroupConfig.T(),
default=[WaveformTargetGroupConfig()])
gnss_config = GNSSTargetGroupConfig.T(
default=GNSSTargetGroupConfig(), optional=True)
insar_config = InSARTargetGroupConfig.T(
default=InSARTargetGroupConfig(), optional=True)
run_config = RunConfig.T(default=RunConfig())
engine_config = EngineConfig.T(default=EngineConfig())
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
gnss_config = kwargs.get('gnss_config')
if gnss_config is not None:
self.gnss_config = GNSSTargetGroupConfigFull(**gnss_config.items())
self.bind_parent()
def bind_parent(self):
for conf in (
self.dataset_config,
self.gnss_config,
self.insar_config,
self.run_config):
if conf is None:
continue
conf.bind_parent(self)
for conf in self.waveform_config:
conf.bind_parent(self)
class IDSConfigFull(IDSConfig):
@property
def inputdir(self):
return 'input'
@property
def inputfn(self):
return op.join(self.inputdir, 'ids_input.inp')
@property
def waveformdir(self):
return op.join(self.inputdir, 'waveform')
@property
def gnssdir(self):
return op.join(self.inputdir, 'gnss')
@property
def gnssfile(self):
fn = 'gnss_data.dat'
if self.gnss_config:
fn = self.gnss_config.fn_dat
return op.join(self.gnssdir, fn)
@property
def insardir(self):
return op.join(self.inputdir, 'insar')
@property
def rundir(self):
return 'run'
@staticmethod
def example():
dataset_config = DatasetConfig.example()
waveform_config = [WaveformTargetGroupConfig.example()]
gnss_config = GNSSTargetGroupConfig.example()
insar_config = InSARTargetGroupConfig.example()
engine_config = EngineConfig.example()
run_config = RunConfig.example()
return IDSConfigFull(
path_prefix='..',
dataset_config=dataset_config,
waveform_config=waveform_config,
gnss_config=gnss_config,
insar_config=insar_config,
engine_config=engine_config,
run_config=run_config)
def string_for_config(self, dataset):
event = dataset.event
eng_conf = self.engine_config
np = self.dataset_config.fault_plane[-1]
gnss_target = dataset.gnss_targetgroup
d = dict(
origin_time=putil.tts(event.time, format='%Y %m %d %H %M %S'),
lat=event.lat,
lon=event.lon,
depth=event.depth / km,
mag=event.magnitude,
waveform_dir=op.relpath(
self.waveformdir, start=self.inputdir) + os.sep,
dist_min=dataset.get_waveform_distance_range()[0] / km,
dist_max=dataset.get_waveform_distance_range()[1] / km,
waveform_static_weight=dataset.get_waveform_static_weight(),
gf_dir=eng_conf.ids_waveform_store,
filter_order=eng_conf.ids_synthetic_bp_filter.order,
filter_f_min=eng_conf.ids_synthetic_bp_filter.f_min,
filter_f_max=eng_conf.ids_synthetic_bp_filter.f_max,
gnss_nsta=gnss_target.n_stations if gnss_target else 0,
gnss_dat=op.relpath(
self.gnssfile, start=self.inputdir), # TODO
gnss_weight=dataset.gnss_targetgroup.config.weight, # TODO
insar_ngrid=self.insar_config.n_grids, # TODO
insar_dat=self.insar_config.dat_path, # TODO
insar_weight=self.insar_config.weight, # TODO
geodetic_gf_dir=eng_conf.ids_geodetic_store,
geodetic_gf_z=eng_conf.ids_geodetic_files['Z'],
geodetic_gf_radial=eng_conf.ids_geodetic_files['R'],
geodetic_gf_transvers=eng_conf.ids_geodetic_files['T'],
idisc=1,
subfault='#',
subfault_file="' '",
n_subfaults=0,
ref_subfault=None,
focmec='',
strike=getattr(event.moment_tensor, 'strike{}'.format(np)),
dip=getattr(event.moment_tensor, 'dip{}'.format(np)),
rake=getattr(event.moment_tensor, 'rake{}'.format(np)),
n_iter_max=self.run_config.n_iterations,
run_dir=op.relpath(self.rundir, start=self.inputdir) + os.sep)
t_max = self.run_config.t_max_snapshots
deltat = self.run_config.deltat_snapshots
snaps = ["'Snapshot_{}s.dat' {:<10.1f}{:<10.1f}".format(t, 0.0, t)
for t in num.arange(deltat, t_max + deltat, deltat)]
snaps += [
"'SnapshotR_{}s.dat' {:<10.1f}{:<10.1f}".format(t, t - deltat, t)
for t in num.arange(deltat, t_max + deltat, deltat)]
d['n_snaps'] = len(snaps)
d['snapshot_lines'] = '\n '.join(snaps)
template = '''# Example for the input file of FORTRAN77 program "ids2020" for kinematic earthquake
# source imaging based on the Iterative Deconvolution and Stacking (IDS) Method using
# seismic waveform and static displacement data jointly.
#
# written by Rongjiang Wang
# Helmholtz Centre Potsdam
# GFZ German Research Centre for Geosciences
# e-mail: wang@gfz-potsdam.de
# Phone +49 331 2881209
# Fax +49 331 2881204
#
# Last modified: Potsdam, August, 2020
#
# Reference:
# Zhang and Wang et al. (2014). Automatic imaging of earthquake rupture processes
# by iterative deconvolution and stacking of high-rate GNSS and strong motion
# seismograms, JGR Solid Earth, 199(7), 5633-5650.
#
#################################################################
## ##
## If not specified otherwise, SI Unit System is used overall! ##
## ##
## ##
## Each comment line should start with "#" ##
## and no blank line is allowed! ##
## ##
#################################################################
#
#================================================================================
# EARTHQUAKE LOCATION PARAMETERS
# ------------------------------
# 1. origin time [s] (year, month, day, hour, minute [positive integer number]
# and second [positive or negative real number]), which will be also used as
# the reference time of the output data
# 2. hypocentre (lat[deg], lon[deg], depth[km])
# 3. preliminary estimate of moment magnitude
# Note: the hypocentre location and origin time (at least to minute) should be
# consistent with those given in the "StationInfo.dat" file of the
# waveform data (s. below)
#--------------------------------------------------------------------------------
{origin_time}
{lat} {lon} {depth}
{mag}
#================================================================================
# WAVEFORM DATA (HIGH-RATE-GNSS/STRONG-MOTION SEISMOGRAMS)
# --------------------------------------------------------------------
# 1. folder name of the data
# Note: this folder should include a file 'StationInfo.dat' giving all
# necessary information about the local monitoring network.
# 2. lower and upper thresholds of the epicentral distances [km] (stations
# outside the threshold distances will not be used)
# 3. weight of the static offset data (obtained automatically from the velocity
# waveforms after an empirical baseline correction) relative to the waveform
# data (suggested value 0.0-1.0)
#--------------------------------------------------------------------------------
'{waveform_dir}'
{dist_min} {dist_max}
{waveform_static_weight}
#================================================================================
# GREEN'S FUNCTION DATABASE
# -------------------------
# 1. Green's function database (folder name)
# 2. Parameters of Butterworth bandpass filter applied to the synthetics:
# order, lower and upper cutoff frequencies [Hz]
# computation [Hz]
#--------------------------------------------------------------------------------
'{gf_dir}'
{filter_order} {filter_f_min} {filter_f_max}
#================================================================================
# GNSS STATIC OFFSETS
# -------------------
# 1. number of GNSS sites, file name of the data
# Note: if the number of GNSS sites <= 0, then no GNSS data are available.
# the GNSS data file includes one header line and 8 culomns,
# Lat(deg) Lon(deg) East(m) North(m) Up(m) wf_E wf_N wf_Z
# where wf_E/N/Z are within-dataset relative weighting factor proportional to
# inverse error variance
# 2. weight of the whole GNSS dataset relative to the waveform dataset (suggested
# value 0.0-1.0)
#--------------------------------------------------------------------------------
{gnss_nsta} '{gnss_dat}'
{gnss_weight}
#================================================================================
# InSAR LOS DISPLACEMENTS
# -----------------------
# 1. number of InSAR grids, file name of the data
# Note: if the number of InSAR grids <= 0, then no InSAR data are used.
# the InSAR data file includes one header line and 8 culomns,
# Lat(deg) Lon(deg) LOS(m) wf_los incidence(deg) azimuth(deg)
# where wf_los is the within-dataset relative weighting factor
# proportional to inverse error variance
# 2. weight of the whole InSAR dataset in the misfit variance relative to the
# waveform dataset (suggested value 0.0-1.0)
#--------------------------------------------------------------------------------
{insar_ngrid} '{insar_dat}'
{insar_weight}
#================================================================================
# Geodetic GREEN'S FUNCTION DATABASE
# ----------------------------------
# 1. geodetic Green's function database (folder name)
# 2. 3 file names of the geodetic Green's functions
#--------------------------------------------------------------------------------
'{geodetic_gf_dir}'
'{geodetic_gf_z}' '{geodetic_gf_radial}' '{geodetic_gf_transvers}'
#================================================================================
# FAULT PARAMETERS
# ----------------
# 1. ipatch: switch (0/1) for subfaults input form,
# 0 = read discrete patches from files,
# 1 = automatic discretization of a single rectangular fault
#--------------------------------------------------------------------------------
# idisc
#--------------------------------------------------------------------------------
{idisc}
#--------------------------------------------------------------------------------
# IF(idisc = 0)THEN
#
# subfault paramters
# ==================
# file_name: data file for discrete patches
# number of sub-faults
# iref selection of the patch reference location
#
# Note: the file includes one header line followed by 8 columns of data:
#
# latitude[deg], longitude[deg], depth[km], length[km], width[km],
# strike[deg], dip[deg], rake[deg]
#
# where (lat, lon, dep) is the reference location of the patch depending
# on iref: 1 = upper-left corner, 2 = upper-right corner,
# 3 = lower-left corner, 4 = lower-right corner,
# 5 = central point (s. Fig.)
#
# definitions for a rectangular fault patch
# =========================================
#
# north(x)
# /
# / | strike
# 1-----------------------2
# |\ p . \ W
# :-\ i . \ i
# | \ l . 5 \ d
# :90 \ S . \ t
# |-dip\ . \ h
# : \. ) rake \
# downward(z) 3-----------------------4
# L e n g t h
#--------------------------------------------------------------------------------
# file_name nsf iref
#--------------------------------------------------------------------------------
{subfault} {subfault_file} {n_subfaults} {ref_subfault}
#--------------------------------------------------------------------------------
# ELSE IF(idisc = 1)THEN
#
# 2. focal mechanism: strike, dip and rake angles [deg]
#--------------------------------------------------------------------------------
{focmec} {strike:8.3f} {dip:7.3f} {rake:8.3f}
#================================================================================
# IDS PARAMETERS
# --------------
# 1. maximum numbers of iterations.
# Note: if max. iteration = 0, forward modelling for any dummy data set is
# made using the given rupture model and the filter parameters
# if max. iteration < 0, same as iteration = 0, but without filtering.
#--------------------------------------------------------------------------------
{n_iter_max}
#================================================================================
# OUTPUT
# ------
# 1. folder name for the output files (should exist already)
# 2. sub-folder name for comparisons of observed and modelled waveform data
# (should exist already)
# 3. sub-folder name for results from empirical baseline correction of
# waveform data (this folder should exist already)
# 4. file name of major earthquake parameters
# 5. file name of earthquake source-time-function (moment rate history)
# 6. file name of sub-fault source-time-function (moment rate history)
# 7. file name of rupture propagation
# 8. file name of final slip model
# Note: if forward modelling is selected (s. above), the above two files
# should provide data for the given kinematic rupture model and will
# be overwritten after the IDS inversion
# 9. file name of normalized residual waveform variance
#10. file name of predicted static offsets at seismic waveform data sites
#11. file name of comparison between observed and modelled GNSS data
#12. file name of comparison between observed and modelled InSAR data
#13. number of snapshots of the rupture process
#14. file name of the 1. snapshot for the slip released within the given
# time window (t1[s], t2[s])
#--------------------------------------------------------------------------------
'{run_dir}'
'O2S'
'EBC'
'EQ_MP.dat'
'EQ_STF.dat'
'PT_STF.dat'
'Rup_Front.dat'
'SlipModel.dat'
'WFMisfitvar.dat'
'SMStaticOffset.dat'
'GNSSDatafit.dat'
'InSARDatafit.dat'
{n_snaps}
{snapshot_lines}
#================================end of input====================================
''' # noqa
return template.format(**d).encode('ascii')
__all__ = [
'IDSConfig',
'IDSConfigFull',
'DatasetConfig',
'EngineConfig',
'RunConfig']

503
src/si/ids/core.py

<
@ -1,255 +1,248 @@
# GPLv3
#
# The Developers, 21st Century
import logging
import tempfile
import shutil
import signal
from subprocess import Popen, PIPE
import os
import os.path as op
from pyrocko.guts import Object, Bool
from pyrocko.io import save
from ewrica.sources import dump_sources
from ewrica.io.ids import load_ids_source, load_obssyn
from .config import IDSConfigFull
from .dataset import Dataset
from .targets import WaveformTargetGroup, GNSSTargetGroup, InSARTargetGroup
logger = logging.getLogger('ewrica.si.ids.core')
km = 1e3
# unit2num = dict([
# ('M', 0),
# ('M/S', 1),
# ('M/S**2', 2)])
program_bins = {
'ids.2020': 'ids2020'
}
class IDSError(Exception):
pass
class IDSRunner(Object):
config = IDSConfigFull.T(
default=IDSConfigFull())
keep_tmp = Bool.T(default=False)
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self._tempdir = ''
self._dataset = None
# self._event = None
# self._p_arrivals = []
# self._sx = None
# self._stations = []
# self._unique_station_targets = False
@property
def tempdir(self):
if not self._tempdir:
self._tempdir = tempfile.mkdtemp('', 'idsrunner-tmp-')
return self._tempdir
@property
def inputdir(self):
return op.join(self.tempdir, self.config.inputdir)
@property
def waveformdir(self):
return op.join(self.tempdir, self.config.waveformdir)
@property
def gnssdir(self):
return op.join(self.tempdir, self.config.gnssdir)
@property
def gnssfile(self):
return op.join(self.tempdir, self.config.gnssfile)
@property
def insardir(self):
return op.join(self.tempdir, self.config.insardir)
@property
def rundir(self):
return op.join(self.tempdir, self.config.rundir)
def get_dataset(self):
if self._dataset is None:
wf_targets = []
gnss_target = None
insar_target = None
for c in self.config.waveform_config:
wf_targets.append(WaveformTargetGroup(config=c))
if self.config.gnss_config:
gnss_target = GNSSTargetGroup(config=self.config.gnss_config)
if self.config.insar_config:
insar_target = InSARTargetGroup(
config=self.config.insar_config)
self._dataset = Dataset(
config=self.config.dataset_config,
waveform_targetgroups=wf_targets,
gnss_targetgroup=gnss_target,
insar_targetgroup=insar_target)
return self._dataset
def _prepare_ids_dir(self):
dirs = \
[self.waveformdir, self.gnssdir, self.insardir] + \
[op.join(self.rundir, d) for d in ('O2S', 'EBC')]
for d in dirs:
os.makedirs(d)
def prepare(self, *args, **kwargs):
self._prepare_ids_dir(*args, **kwargs)
ds = self.get_dataset()
ds.set_waveform_data(self.waveformdir)
ds.set_gnss_data(self.gnssdir)
ds.set_insar_data(self.insardir)
self.config.dump(filename=op.join(
self.tempdir, self.config.inputdir, 'ids_config.yaml'))
def run(self):
config = self.config
input_fn = op.join(self.tempdir, config.inputfn)
self._input_fn = input_fn
with open(input_fn, 'wb') as f:
input_str = config.string_for_config(self.get_dataset())
logger.debug('===== begin ids input =====\n'
'%s===== end ids input =====' % input_str.decode())
f.write(input_str)
program = program_bins['ids.%s' % config.ids_version]
old_wd = os.getcwd()
os.chdir(self.inputdir)
interrupted = []
def signal_handler(signum, frame):
os.kill(proc.pid, signal.SIGTERM)
interrupted.append(True)
original = signal.signal(signal.SIGINT, signal_handler)
try:
try:
proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
except OSError:
os.chdir(old_wd)
raise IDSError(
'''could not start ids executable: "%s"
Available ids backends and download links to the modelling codes are listed
on
https://git.pyrocko.org/wanderelch/ewrica
''' % program)
(output_str, error_str) = proc.communicate(
'{}\n'.format(input_fn).encode())
finally:
signal.signal(signal.SIGINT, original)
if interrupted:
raise KeyboardInterrupt()
logger.info('===== begin ids output =====\n'
'%s===== end ids output =====' % output_str.decode())
errmess = []
if proc.returncode != 0:
errmess.append(
'ids had a non-zero exit state: %i' % proc.returncode)
if error_str:
if error_str.lower().find(b'error') != -1:
errmess.append("the string 'error' appeared in ids error")
else:
logger.warn(
'ids emitted something via stderr:\n\n%s'
% error_str.decode())
if output_str.lower().find(b'error') != -1:
errmess.append("the string 'error' appeared in ids output")
if errmess:
self.keep_tmp = True
os.chdir(old_wd)
raise IDSError('''
===== begin ids input =====
%s===== end ids input =====
===== begin ids output =====
%s===== end ids output =====
===== begin ids error =====
%s===== end ids error =====
%s
ids has been invoked as "%s"
in the directory %s'''.lstrip() % (
input_str.decode(), output_str.decode(), error_str.decode(),
'\n'.join(errmess), program, self.tempdir))
self.ids_output = output_str
self.ids_error = error_str
os.chdir(old_wd)
def post_process(self, force=False):
source = load_ids_source(self._input_fn)
outpath_tmpl = self.config.run_config.outpath_template
if outpath_tmpl is not None:
result_dir = op.join(self.tempdir, 'results')
os.makedirs(result_dir)
dump_sources(
[source],
filename=op.join(result_dir, 'idssource.yaml'))
traces = load_obssyn(
config_fn=op.join(self.tempdir, self.config.inputfn))
save(
traces, filename_template=op.join(result_dir, 'traces.mseed'))
try:
shutil.copytree(
self.tempdir,
outpath_tmpl,
dirs_exist_ok=force)
except OSError:
raise IDSError('Configured outpath is already used. Consider '
'using the "force" keyword argument')
if not self.keep_tmp:
shutil.rmtree(self.tempdir)
return source
# GPLv3
#
# The Developers, 21st Century
import logging
import tempfile
import shutil
import signal
from subprocess import Popen, PIPE
import os
import os.path as op
from pyrocko.guts import Object, Bool
from pyrocko.io import save
from ewrica.sources import dump_sources
from ewrica.io.ids import load_ids_source, load_obssyn
from .config import IDSConfigFull
from .dataset import Dataset
from .targets import WaveformTargetGroup, GNSSTargetGroup, InSARTargetGroup
logger = logging.getLogger('ewrica.si.ids.core')
km = 1e3
program_bins = {
'ids.2020': 'ids2020'
}
class IDSError(Exception):
pass
class IDSRunner(Object):
config = IDSConfigFull.T(
default=IDSConfigFull())
keep_tmp = Bool.T(default=False)
def __init__(self, *args, **kwargs):
Object.__init__(self, *args, **kwargs)
self._tempdir = ''
self._dataset = None
@property
def tempdir(self):
if not self._tempdir:
self._tempdir = tempfile.mkdtemp('', 'idsrunner-tmp-')
return self._tempdir
@property
def inputdir(self):
return op.join(self.tempdir, self.config.inputdir)
@property
def waveformdir(self):
return op.join(self.tempdir, self.config.waveformdir)
@property
def gnssdir(self):
return op.join(self.tempdir, self.config.gnssdir)
@property
def gnssfile(self):
return op.join(self.tempdir, self.config.gnssfile)
@property
def insardir(self):
return op.join(self.tempdir, self.config.insardir)
@property
def rundir(self):
return op.join(self.tempdir, self.config.rundir)
def get_dataset(self):
if self._dataset is None:
wf_targets = []
gnss_target = None
insar_target = None
for c in self.config.waveform_config:
wf_targets.append(WaveformTargetGroup(config=c))
if self.config.gnss_config:
gnss_target = GNSSTargetGroup(config=self.config.gnss_config)
if self.config.insar_config:
insar_target = InSARTargetGroup(
config=self.config.insar_config)
self._dataset = Dataset(
config=self.config.dataset_config,
waveform_targetgroups=wf_targets,
gnss_targetgroup=gnss_target,
insar_targetgroup=insar_target)
return self._dataset
def _prepare_ids_dir(self):
dirs = \
[self.waveformdir, self.gnssdir, self.insardir] + \
[op.join(self.rundir, d) for d in ('O2S', 'EBC')]
for d in dirs:
os.makedirs(d)
def prepare(self, *args, **kwargs):
self._prepare_ids_dir(*args, **kwargs)
ds = self.get_dataset()
ds.set_waveform_data(self.waveformdir)
ds.set_gnss_data(self.gnssdir)
ds.set_insar_data(self.insardir)
self.config.dump(filename=op.join(
self.tempdir, self.config.inputdir, 'ids_config.yaml'))
def run(self):
config = self.config
input_fn = op.join(self.tempdir, config.inputfn)
self._input_fn = input_fn
with open(input_fn, 'wb') as f:
input_str = config.string_for_config(self.get_dataset())
logger.debug('===== begin ids input =====\n'
'%s===== end ids input =====' % input_str.decode())
f.write(input_str)
program = program_bins['ids.%s' % config.ids_version]
old_wd = os.getcwd()
os.chdir(self.inputdir)
interrupted = []
def signal_handler(signum, frame):
os.kill(proc.pid, signal.SIGTERM)
interrupted.append(True)
logger.info('Start {} computation'.format(program))
original = signal.signal(signal.SIGINT, signal_handler)
try:
try:
proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
except OSError:
os.chdir(old_wd)
raise IDSError(
'''could not start ids executable: "%s"
Available ids backends and download links to the modelling codes are listed
on
https://git.pyrocko.org/wanderelch/ewrica
''' % program)
(output_str, error_str) = proc.communicate(
'{}\n'.format(input_fn).encode())
finally:
signal.signal(signal.SIGINT, original)
if interrupted:
raise KeyboardInterrupt()
logger.info('===== begin ids output =====\n'
'%s===== end ids output =====' % output_str.decode())
errmess = []
if proc.returncode != 0:
errmess.append(
'ids had a non-zero exit state: %i' % proc.returncode)
if error_str:
if error_str.lower().find(b'error') != -1:
errmess.append("the string 'error' appeared in ids error")
else:
logger.warn(
'ids emitted something via stderr:\n\n%s'
% error_str.decode())
if output_str.lower().find(b'error') != -1:
errmess.append("the string 'error' appeared in ids output")
if errmess:
self.keep_tmp = True
os.chdir(old_wd)
raise IDSError('''
===== begin ids input =====
%s===== end ids input =====
===== begin ids output =====
%s===== end ids output =====
===== begin ids error =====
%s===== end ids error =====
%s
ids has been invoked as "%s"
in the directory %s'''.lstrip() % (
input_str.decode(), output_str.decode(), error_str.decode(),
'\n'.join(errmess), program, self.tempdir))
self.ids_output = output_str
self.ids_error = error_str
os.chdir(old_wd)
def post_process(self, force=False):
source = load_ids_source(self._input_fn)
outpath_tmpl = self.config.run_config.outpath_template
if outpath_tmpl is not None:
result_dir = op.join(self.tempdir, 'results')
os.makedirs(result_dir)
dump_sources(
[source],
filename=op.join(result_dir, 'idssource.yaml'))
traces = load_obssyn(
config_fn=op.join(self.tempdir, self.config.inputfn))
save(
traces, filename_template=op.join(result_dir, 'traces.mseed'))
if op.exists(outpath_tmpl) and not force:
raise IDSError('Configured outpath is already used. Consider '
'using the "force" keyword argument')
if op.exists(outpath_tmpl):
shutil.rmtree(outpath_tmpl)
shutil.copytree(
self.tempdir,
outpath_tmpl,
dirs_exist_ok=force)