Browse Source

ids: wrapper for fortran IDS function

- full wrapping for
	- waveform loading
	- restitution
	- DAT file saving
	- ids fortran code call
	- result and output handling
- WIP:
	- tests
	- loading of synthetic DAT files for possible comparison
mole
mmetz 10 months ago
parent
commit
8b72325977
  1. 4
      setup.py
  2. 989
      src/io/ids.py
  3. 0
      src/si/ids/__init__.py
  4. 502
      src/si/ids/config.py
  5. 612
      src/si/ids/core.py
  6. 7
      src/si/ids/targets/__init__.py
  7. 78
      src/si/ids/targets/base.py
  8. 21
      src/si/ids/targets/gnss.py
  9. 21
      src/si/ids/targets/insar.py
  10. 216
      src/si/ids/targets/waveform.py
  11. 16
      src/sources.py
  12. 20
      src/util.py
  13. 0
      test/io/__init__.py
  14. 712
      test/io/test_ids.py
  15. 0
      test/si/__init__.py
  16. 63
      test/si/test_ids.py

4
setup.py

@ -49,7 +49,9 @@ packages = [
'%s.ls' % packname,
'%s.message' % packname,
'%s.plot' % packname,
'%s.si' % packname]
'%s.si' % packname,
'%s.si.ids' % packname,
'%s.si.ids.targets' % packname]
entry_points = {

989
src/io/ids.py

File diff suppressed because it is too large

0
src/data/examples → src/si/ids/__init__.py

502
src/si/ids/config.py

@ -0,0 +1,502 @@
# GPLv3
#
# The Developers, 21st Century
import logging
import os
import os.path as op
import numpy as num
from pyrocko import model, util as putil
from pyrocko.guts import Float, List, String, Dict, StringChoice, Int, Object
from .targets import \
TargetConfig, WaveformTargetConfig, GNSSTargetConfig, InSARTargetConfig
logger = logging.getLogger('ewrica.si.ids.config')
km = 1e3
class DatasetConfig(Object):
events_path__ = String.T(default='event.txt')
event_name = String.T(default='event1')
fault_plane = StringChoice.T(
choices=['nodalplane1', 'nodalplane2'], default='nodalplane1')
_expand_paths = TargetConfig._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 DatasetConfig()
@property
def events_path(self):
return self._expand_paths(self.events_path__)
@events_path.setter
def events_path(self, events_path):
self.events_path__ = events_path
@property
def event(self):
events = model.load_events(self.events_path)
for event in events:
if event.name == self.event_name:
return event
break
raise ValueError(
'Given event_name "{}" not found in given event file'.format(
self.event_name))
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 = TargetConfig._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_conf = DatasetConfig.T(default=DatasetConfig())
waveform_conf = WaveformTargetConfig.T(
default=WaveformTargetConfig())
gnss_conf = GNSSTargetConfig.T(
default=GNSSTargetConfig(), optional=True)
insar_conf = InSARTargetConfig.T(
default=InSARTargetConfig(), optional=True)
run_conf = RunConfig.T(default=RunConfig())
engine_conf = 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_conf,
self.waveform_conf,
self.gnss_conf,
self.insar_conf,
self.run_conf):
if conf is None:
continue
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 insardir(self):
return op.join(self.inputdir, 'insar')
@property
def rundir(self):
return 'run'
@staticmethod
def example():
dataset_conf = DatasetConfig.example()
waveform_conf = WaveformTargetConfig.example()
gnss_conf = GNSSTargetConfig.example()
insar_conf = InSARTargetConfig.example()
engine_conf = EngineConfig.example()
run_conf = RunConfig.example()
return IDSConfigFull(
path_prefix='..',
dataset_conf=dataset_conf,
waveform_conf=waveform_conf,
gnss_conf=gnss_conf,
insar_conf=insar_conf,
engine_conf=engine_conf,
run_conf=run_conf)
def string_for_config(self):
event = self.dataset_conf.event
eng_conf = self.engine_conf
np = self.dataset_conf.fault_plane[-1]
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=self.waveform_conf.distance_min / km,
dist_max=self.waveform_conf.distance_max / km,
waveform_static_weight=self.waveform_conf.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=self.gnss_conf.n_stations, # TODO
gnss_dat=self.gnss_conf.dat_path, # TODO
gnss_weight=self.gnss_conf.weight, # TODO
insar_ngrid=self.insar_conf.n_grids, # TODO
insar_dat=self.insar_conf.dat_path, # TODO
insar_weight=self.insar_conf.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_conf.n_iterations,
run_dir=op.relpath(self.rundir, start=self.inputdir) + os.sep)
t_max = self.run_conf.t_max_snapshots
deltat = self.run_conf.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} {dip} {rake}
#================================================================================
# 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']

612
src/si/ids/core.py

@ -0,0 +1,612 @@
# 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
import numpy as num
from pyrocko import marker, gf, util as putil, pile, orthodrome as pod, model
from pyrocko.guts import Object, Bool
from pyrocko.io import stationxml
from ewrica import util
from ewrica.sources import dump_sources
from ewrica.io.ids import load_ids_source
from .config import IDSConfigFull
from .targets import RestitutionConfig
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 RestitutionTarget(Object):
station = model.Station.T(optional=True)
restitution_conf = RestitutionConfig.T(optional=True)
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._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 insardir(self):
return op.join(self.tempdir, self.config.insardir)
@property
def rundir(self):
return op.join(self.tempdir, self.config.rundir)
@property
def event(self):
if self._event is None:
self._event = self.config.dataset_conf.event
return self._event
@property
def sx(self):
wf_conf = self.config.waveform_conf
if self._sx is None:
sxs = []
for f in wf_conf.response_stationxml_paths:
sxs.append(stationxml.load_xml(filename=f))
self._sx = stationxml.primitive_merge(sxs)
return self._sx
@property
def stations(self):
if not self._stations:
wf_conf = self.config.waveform_conf
stations = self.sx.get_pyrocko_stations()
stations = util.keep_allowlisted_stations(
stations, wf_conf.allowed)
stations = util.remove_blocklisted_stations(
stations, wf_conf.blocked)
self._stations = stations
return self._stations
@property
def event_station_distance(self):
return pod.distance_accurate50m_numpy(
self.event.lat, self.event.lon,
[s.lat for s in self.stations], [s.lon for s in self.stations])
@property
def p_arrivals(self):
if self._p_arrivals:
return self._p_arrivals
wf_conf = self.config.waveform_conf
en_conf = self.config.engine_conf
self._ensure_unique_station_targets()
arrivals = []
for target in wf_conf.waveform_targets:
if (target.ppicks_paths is not None and target.pphase is None and
target.store_id is None):
for fn in target.ppicks_paths:
arrivals += marker.load_markers(fn)
elif (target.ppicks_paths is None and target.pphase is not None and
target.store_id is not None):
engine = gf.LocalEngine(
store_superdirs=en_conf.pyrocko_store_superdirs)
store = engine.get_store(target.store_id)
try:
store_arrs = num.array([
store.t(target.pphase, self.event, sta)
for sta in self.stations]) + self.event.time
except gf.store.NoSuchPhase as e:
raise e
markers = [
marker.Marker(
tmin=arr, tmax=arr,
nslc_ids=[s.nsl() + (c,) for c in ['*E', '*N', '*Z']])
for arr, s in zip(store_arrs, self.stations)]
arrivals += markers
else:
raise ValueError('Argument "ppicks" is mutually exclusive '
'with "phase" and "store_id".')
self._p_arrivals = arrivals
return self._p_arrivals
def _ensure_unique_station_targets(self):
if self._unique_station_targets:
pass
stations = self.stations
targets = self.config.waveform_conf.waveform_targets
stations_used = []
for target in targets:
stations_relevant = [
s for s, d in zip(self.stations, self.event_station_distance)
if d >= target.distance_min and d < target.distance_max]
stations_relevant = util.keep_allowlisted_stations(
stations, target.allowed)
stations_relevant = util.remove_blocklisted_stations(
stations_relevant, target.blocked)
for s in stations_relevant:
if s.nsl() in stations_used:
raise IDSError(
'{} - station traces may only be allowed in one '
'target groups'.format('.'.join(s.nsl())))
stations_used += stations_relevant
target._stations = stations_relevant
self._unique_station_targets = True
def _generate_stationinfo(self, station_tuples):
# TODO update for waveform target list
stations = self.stations
nsls = [st[0] for st in station_tuples]
t_pres = [st[1] for st in station_tuples]
t_ps = [st[2] for st in station_tuples]
deltats = [st[3] for st in station_tuples]
n_samples = [st[4] for st in station_tuples]
units = [st[5] for st in station_tuples]
weights = [st[6] for st in station_tuples]
t_base = self.event.time - self.event.time % 60
t_min = num.min(t_pres + [self.event.time])
order = 'E N Z'.split()
order2num = dict(E=1, N=2, Z=3)
template = '''# reference time: year month day hour minute second
# (year - minute in integer format with the usual convention, second by any real value,
# e.g., 0.0, 25.0, -30.0, 75.0, ...)
{reference_time}
# earthquake hypocentre: lat[deg] lon[deg] depth[km]
{lat} {lon} {depth}
# number of stations
{n_stations}
# column of component
# east north up => 1 2 3
{col_1:d} {col_2:d} {col_3:d}
# Note: Start and P_onset are relative to the reference time, T_length is the usable length of data,
# Weight_1-3 are the relative weighting factors of the 3 components of data
# P_onset <= 0: not available
# T_length[s] <= 0: not available
# Data_type: 0/1/2 = displacement/velocity/acceleration
# Sampling_interval: station dependent sampling interval
#
# Station Lat[deg] Lon[deg] Start[s] P_onset[s] T_length[s] Weight_1 Weight_2 Weight_3 Data_type Sampling_interval[s]
{station_lines:s}
''' # noqa
d = dict(
reference_time='%s %.1f' % (
putil.tts(t_base, format='%Y %m %d %H %M'), t_min - t_base),
lat=self.event.lat,
lon=self.event.lon,
depth=self.event.depth / km,
col_1=order2num[order[0]],
col_2=order2num[order[1]],
col_3=order2num[order[2]])
lines = []
nsta = 0
line_fmt = (
'{station:<15}{lat:<15g}{lon:<15g}{tmin:<15g}{tp:<15g}'
'{tlength:<15g}{weight1:<15g}{weight2:<15g}{weight3:<15g}'
'{unit:<15}{deltat:<15g}')
stations_used = {}
for sta in stations:
nsl = '.'.join([x for x in sta.nsl() if x])
if nsl not in nsls:
continue
if nsl in stations_used:
logger.warn(
'{} - Multiple station information. Used \n{}\n\n'.format(
nsl, sta.__str__()))
idx = nsls.index(nsl)
lines.append(line_fmt.format(
station=sta.station,
lat=sta.lat,
lon=sta.lon,
tmin=t_pres[idx] - t_min,
tp=t_ps[idx] - t_min,
tlength=n_samples[idx] * deltats[idx],
weight1=weights[idx][order[0]],
weight2=weights[idx][order[1]],
weight3=weights[idx][order[2]],
unit=unit2num[units[idx]],
deltat=deltats[idx]))
nsta += 1
stations_used[nsl] = sta
d['station_lines'] = '\n '.join(lines)
d['n_stations'] = nsta
with open(op.join(self.waveformdir, 'StationInfo.dat'), 'wb') as f:
f.write((template.format(**d)).encode('ascii'))
def _restitute_waveform_data(self, *args, **kwargs):
# TODO update for waveform target list
config = self.config
wf_conf = config.waveform_conf
p = pile.make_pile(wf_conf.waveform_paths)
traces_rest = []
fake_input_units = []
weights = []
for target in wf_conf.waveform_targets:
rest_conf = target.restitution_conf
time_range = (
self.event.time + rest_conf.t_min,
self.event.time + rest_conf.t_max)
t_fade = rest_conf.t_fade
if t_fade is None:
t_fade = 1. / num.min(rest_conf.frequency_limits)
stations_relevant = {s.nsl(): s for s in target._stations}
for traces in p.chopper_grouped(
gather=lambda tr: tr.nslc_id, want_incomplete=False,
tmin=time_range[0]-t_fade, tmax=time_range[1]+t_fade):
for tr in traces:
if tr.nslc_id[:3] not in stations_relevant:
continue
try:
polezero_response = self.sx.get_pyrocko_response(
nslc=tr.nslc_id,
timespan=(tr.tmin, tr.tmax),
fake_input_units=rest_conf.target_unit)
except stationxml.NoResponseInformation:
logger.warn(
'{} excluded - No reponse found'.format(
'.'.join(tr.nslc_id)))
continue
except stationxml.MultipleResponseInformation:
logger.warn(
'{} excluded - multiple reponses found'.format(
'.'.join(tr.nslc_id)))
continue
traces_rest.append(tr.transfer(
tfade=t_fade,
freqlimits=rest_conf.frequency_limits,
transfer_function=polezero_response,
invert=True).chop(
tmin=time_range[0], tmax=time_range[1]))
if rest_conf.target_unit is None:
unit_resp = [
c.response.instrument_sensitivity.input_units.name
for n, s, c in
self.sx.iter_network_station_channels(
net=tr.nslc_id[0],
sta=tr.nslc_id[1],
loc=tr.nslc_id[2])]
if num.unique(unit_resp).shape[0] > 1:
raise IDSError(
'{} - multiple input units found'.format(
'.'.join(tr.nslc_id[:3])))
fake_input_units.append(unit_resp[0])
else:
fake_input_units.append(rest_conf.target_unit)
weights.append(target.component_weights)
if logger.level == logging.DEBUG:
from pyrocko import trace
trace.snuffle(
traces_rest,
stations=self.stations,
markers=self.p_arrivals,
events=[self.event])
return traces_rest, fake_input_units, weights
def _prepare_waveform_data(self, *args, **kwargs):
# TODO update for waveform target list
config = self.config
wf_conf = config.waveform_conf
order = ['E', 'N', 'Z']
self._ensure_unique_station_targets()
traces_rest, input_units, weights = self._restitute_waveform_data()
station_traces = {}
station_units = {}
station_weights = {}
for i, tr in enumerate(traces_rest):
nsl = '.'.join([x for x in tr.nslc_id[:3] if x])
if nsl not in station_traces:
station_traces[nsl] = []
station_traces[nsl].append(tr)
station_units[nsl] = input_units[i]
station_weights[nsl] = weights[i]
arrivals = self.p_arrivals
arrivals = {'.'.join([x for x in ar.nslc_ids[0][:3] if x]): ar
for ar in arrivals}
stations_dat = []
for nsl_string, traces in station_traces.items():
logger.info('{} - Convert to dat'.format(nsl_string))
try:
arrival = arrivals[nsl_string]
except KeyError:
logger.warn('{} - No P arrival times found'.format(nsl_string))
continue
dt = num.unique([t.deltat for t in traces])
if dt.shape[0] > 1:
logger.warn('{} - Different sampling rates'.format(','.join(
['.'.join(t.nslc_id) for t in traces])))
continue
dt = dt[0]
t_p = num.mean([arrival.tmin, arrival.tmax])
t_pre = t_p + wf_conf.t_prearrival
t_pre -= t_pre % dt
tmi = t_pre
tma = num.unique([t.tmax for t in traces])
if len(tma) > 1:
logger.warn('{} - Different trace endtimes'.format(','.join(
['.'.join(t.nslc_id) for t in traces])))
tma = [tma.min()]
tma = tma[0]
traces = [t.chop(tmin=tmi, tmax=tma) for t in traces]
data = num.zeros((traces[0].ydata.shape[0], len(order)))
idx = [order.index(tr.nslc_id[3][-1]) for tr in traces]
for itr, i in enumerate(idx):
data[:, i] = traces[itr].ydata
fn = op.join(self.waveformdir, nsl_string.split('.')[1] + '.dat')
num.savetxt(fn, data)
stations_dat.append(
(nsl_string,
t_pre,
t_p,
dt,
data.shape[0],
station_units[nsl],
station_weights[nsl]))
self._generate_stationinfo(stations_dat)
def _prepare_gnss_data(self, *args, **kwargs):
pass
def _prepare_insar_data(self, *args, **kwargs):
pass
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)
self._prepare_waveform_data(*args, **kwargs)
self._prepare_gnss_data(*args, **kwargs)
self._prepare_insar_data(*args, **kwargs)
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()
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_conf.outpath_template
if outpath_tmpl is not None:
dump_sources(
[source],
filename=op.join(self.rundir, 'idssource.yaml'))
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

7
src/si/ids/targets/__init__.py

@ -0,0 +1,7 @@
# GPLv3
#
# The Developers, 21st Century
from .base import * # noqa
from .waveform import * # noqa
from .gnss import * # noqa
from .insar import * # noqa

78
src/si/ids/targets/base.py

@ -0,0 +1,78 @@
# GPLv3
#
# The Developers, 21st Century
import logging
import os.path as op
from pyrocko.guts import Object, List, String
logger = logging.getLogger('ewrica.si.ids.targets.base')
class TargetError(Exception):
pass
class TargetConfig(Object):
blocklist_paths__ = List.T(String.T(), optional=True)
blocklist = List.T(String.T(), optional=True)
allowlist_paths__ = List.T(String.T(), optional=True)
allowlist = List.T(String.T(), optional=True)
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
def _expand_paths(self, paths):
parent = self._parent
if paths is None:
return None
if parent is None:
return paths
if parent.path_prefix is None or parent.dataset_conf is None:
return paths
event_name = parent.dataset_conf.event_name
if isinstance(paths, list):
return [op.join(
parent.path_prefix,
p.format(event_name=event_name))for p in paths]
elif isinstance(paths, str):
return op.join(
parent.path_prefix, paths.format(event_name=event_name))
else:
raise TargetError(
'{} - Path(s) not formatted as a list or a string'.format(
paths))
@staticmethod
def example():
raise NotImplementedError('Method not implemented')
@property
def blocklist_paths(self):
return self._expand_paths(self.blocklist_paths__)
@blocklist_paths.setter
def blocklist_paths(self, blocklist_paths):
self.blocklist_paths__ = blocklist_paths
@property
def allowlist_paths(self):
return self._expand_paths(self.allowlist_paths__)
@allowlist_paths.setter
def allowlist_paths(self, allowlist_paths):
self.allowlist_paths__ = allowlist_paths

21
src/si/ids/targets/gnss.py

@ -0,0 +1,21 @@
# GPLv3
#
# The Developers, 21st Century
import logging
from pyrocko.guts import Float, List, Int, String
from .base import TargetConfig
logger = logging.getLogger('ewrica.si.ids.targets.gnss')
class GNSSTargetConfig(TargetConfig):
weight = Float.T(default=0.)
gnss_paths = List.T(String.T(), optional=True)
n_stations = Int.T(default=0)
dat_path = String.T(default='./NoGPSData.dat')
@staticmethod
def example():
return GNSSTargetConfig()

21
src/si/ids/targets/insar.py

@ -0,0 +1,21 @@
# GPLv3
#
# The Developers, 21st Century
import logging
from pyrocko.guts import Float, List, Int, String
from .base import TargetConfig
logger = logging.getLogger('ewrica.si.ids.targets.insar')
class InSARTargetConfig(TargetConfig):
weight = Float.T(default=0.)
insar_paths = List.T(String.T(), optional=True)
n_grids = Int.T(default=0)
dat_path = String.T(default='./NoInSARData.dat')
@staticmethod
def example():
return InSARTargetConfig()

216
src/si/ids/targets/waveform.py

@ -0,0 +1,216 @@
# GPLv3
#
# The Developers, 21st Century
import copy
import logging
import numpy as num
from pyrocko.guts import Float, List, String, Tuple, Dict, StringChoice, Object
from .base import TargetConfig
logger = logging.getLogger('ewrica.si.ids.targets.waveform')
km = 1e3
quantity2unit = dict([
('displacement', 'M'),
('velocity', 'M/S'),
('acceleration', 'M/S**2')])
class WaveformQuantities(StringChoice):
choices = ['displacement', 'velocity', 'acceleration']
class RestitutionConfig(Object):
frequency_limits = Tuple.T(4, Float.T(), default=(0.001, 0.05, 100., 200.))
t_fade = Float.T(optional=True)
t_min = Float.T(
default=-60.,
help='start time of raw waveforms relative event origin time')
t_max = Float.T(
default=100.,
help='end time of raw waveforms relative event origin time')
target_quantity = WaveformQuantities.T(optional=True)
@staticmethod
def example():
return RestitutionConfig()
@property
def target_unit(self):
if self.target_quantity is not None:
return quantity2unit[self.target_quantity]
return None
class WaveformTargetGroupConfig(TargetConfig):
distance_min = Float.T(default=0.)
distance_max = Float.T(default=100. * km)
ppicks_paths__ = List.T(String.T(), optional=True)
pphase = String.T(default='{stored:any_P}', optional=True)
store_id = String.T(default='ak135_2000km_1Hz', optional=True)
component_weights = Dict.T(
String.T(), Float.T(),
default=dict(E=0.0, N=0.0, Z=1.0))
restitution_conf = RestitutionConfig.T(default=RestitutionConfig())
def __init__(self, *args, **kwargs):
TargetConfig.__init__(self, *args, **kwargs)
self._blocked = []
self._allowed = []
self._stations = []
@staticmethod
def example():
return WaveformTargetGroupConfig(
blocklist_paths=['path/to/blocklist.txt'],
blocklist=['AB.BLOCK'],
allowlist_paths=['path/to/allowlist.txt'],
allowlist=['AB.ALLOW'],
restitution_conf=RestitutionConfig.example())
@property
def blocked(self):
if not self._blocked:
blocklist = copy.deepcopy(self.blocklist) or []
if self.blocklist_paths:
blocklist += [
str(s) for p in self.blocklist_paths
for s in num.loadtxt(p, dtype='<U20')]
self._blocked = list(num.unique(blocklist))
return self._blocked
@property
def allowed(self):
if not self._allowed:
allowlist = copy.deepcopy(self.allowlist) or []
if self.allowlist_paths:
allowlist += [
str(s) for p in self.allowlist_paths
for s in num.loadtxt(p, dtype='<U20')]
self._allowed = list(num.unique(allowlist))
return self._allowed
@property
def ppicks_path(self):
return self._expand_paths(self.ppicks_path__)
@ppicks_path.setter
def ppicks_path(self, ppicks_path):
self.ppicks_path__ = ppicks_path
class WaveformTargetConfig(TargetConfig):
waveform_paths__ = List.T(String.T(), default=['.'])
response_stationxml_paths__ = List.T(String.T(), default=['response.xml'])
t_prearrival = Float.T(
default=-15.,
help='data starting from time in [sec] relative to P arrivals is '
'included in dat files for baseline corrections.')
static_weight = Float.T(default=0.)
waveform_targets__ = List.T(
WaveformTargetGroupConfig.T(),
default=[WaveformTargetGroupConfig()])
def __init__(self, *args, **kwargs):
TargetConfig.__init__(self, *args, **kwargs)
self._blocked = []
self._allowed = []
@staticmethod
def example():
return WaveformTargetConfig(
waveform_paths=[
'relative/path/to/waveforms1', 'relative/path/to/waveforms1'],
response_stationxml_paths=[
'path/to/responsefile1.xml', 'path/to/responsefile2.xml'],
waveform_targets=[WaveformTargetGroupConfig.example()])
@property
def blocked(self):
if not self._blocked:
blocklist = self.blocklist or []
if self.blocklist_paths:
blocklist += [
str(s) for p in self.blocklist_paths
for s in num.loadtxt(p, dtype='<U20')]
self._blocked = list(num.unique(blocklist))
return self._blocked
@property
def allowed(self):
if not self._allowed:
allowlist = self.allowlist or []
if self.allowlist_paths:
allowlist += [
str(s) for p in self.allowlist_paths
for s in num.loadtxt(p, dtype='<U20')]
self._allowed = list(num.unique(allowlist))
return self._allowed
@property
def waveform_targets(self):
if self._parent is None:
return self.waveform_targets__
if self._parent.path_prefix is None:
return self.waveform_targets__
for t in self.waveform_targets__:
t._parent = self._parent
return self.waveform_targets__
@waveform_targets.setter
def waveform_targets(self, waveform_targets):
self.waveform_targets__ = waveform_targets
@property
def waveform_paths(self):
return self._expand_paths(self.waveform_paths__)
@waveform_paths.setter
def waveform_paths(self, waveform_paths):
self.waveform_paths__ = waveform_paths
@property
def response_stationxml_paths(self):
return self._expand_paths(self.response_stationxml_paths__)
@response_stationxml_paths.setter
def response_stationxml_paths(self, response_stationxml_paths):
self.response_stationxml_paths__ = response_stationxml_paths
@property
def distance_min(self):
return num.min([t.distance_min for t in self.waveform_targets])
@property
def distance_max(self):
return num.min([t.distance_max for t in self.waveform_targets])

16
src/sources.py

@ -296,6 +296,10 @@ class EwricaIDSPatch(gf.Location):
clone = clone
time = Timestamp.T(
optional=True,
help='Rupture front arrival/rupture origin time')
magnitude__ = Float.T(
optional=True,
help='Moment magnitude Mw as in [Hanks and Kanamori, 1979].')
@ -308,10 +312,6 @@ class EwricaIDSPatch(gf.Location):
optional=True,
help='Source time function (moment release) of the patch in [Nm/s].')
time__ = Timestamp.T(
optional=True,
help='Timestamp of the patch rupture start in [s].')
slip = Float.T(
optional=True,
help='Length of the slip vector in [m].')
@ -367,14 +367,6 @@ class EwricaIDSPatch(gf.Location):
if magnitude is not None:
self.magnitude = magnitude
@property
def time(self):
return self._time or num.nan
@time.setter
def time(self, time):
self._time = time
@property
def stf(self):
return self.stf__

20
src/util.py

@ -33,7 +33,10 @@ def setup_file_logging(logger, fn):
logger.addHandler(file_handler)
def remove_blocklisted_stations(stations, stations_blocklist):
def remove_blocklisted_stations(stations, stations_blocklist=[]):
if not stations_blocklist:
return stations
nsl = [s.nsl() for s in stations]
not_blocked = [i for i in range(len(stations))
if nsl[i][1] not in stations_blocklist and
@ -43,7 +46,10 @@ def remove_blocklisted_stations(stations, stations_blocklist):
return [stations[i] for i in not_blocked]
def remove_blocklisted_traces(traces, traces_blocklist):
def remove_blocklisted_traces(traces, traces_blocklist=[]):
if not traces_blocklist:
return traces
nslc = [t.nslc_id for t in traces]
not_blocked = [i for i in range(len(traces))
if nslc[i][1] not in traces_blocklist and
@ -54,7 +60,10 @@ def remove_blocklisted_traces(traces, traces_blocklist):
return [traces[i] for i in not_blocked]
def keep_allowsted_stations(stations, stations_allowlist):
def keep_allowlisted_stations(stations, stations_allowlist=[]):
if not stations_allowlist:
return stations
nsl = [s.nsl() for s in stations]
allowed = [i for i in range(len(stations))
if nsl[i][1] in stations_allowlist or
@ -64,7 +73,10 @@ def keep_allowsted_stations(stations, stations_allowlist):
return [stations[i] for i in allowed]
def keep_allowlisted_traces(traces, traces_allowlist):
def keep_allowlisted_traces(traces, traces_allowlist=[]):
if not traces_allowlist:
return traces
nslc = [t.nslc_id for t in traces]
allowed = [i for i in range(len(traces))
if nslc[i][1] in traces_allowlist or

0
test/io/__init__.py

712
test/io/test_ids.py

@ -0,0 +1,712 @@
import tempfile
import unittest