1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import numpy as num
8import logging
9import os
10import shutil
11import math
12import copy
13import signal
15from tempfile import mkdtemp
16from subprocess import Popen, PIPE
17from os.path import join as pjoin
19from pyrocko import gf
20from pyrocko import trace, util, cake
21from pyrocko.moment_tensor import MomentTensor, symmat6
22from pyrocko.guts import Float, Int, Tuple, List, Complex, Bool, Object, String
24km = 1e3
26guts_prefix = 'pf'
28Timing = gf.meta.Timing
30logger = logging.getLogger('pyrocko.fomosto.qseis')
32# how to call the programs
33program_bins = {
34 'qseis.2006': 'fomosto_qseis2006',
35 'qseis.2006a': 'fomosto_qseis2006a',
36 'qseis.2006b': 'fomosto_qseis2006b',
37}
40def have_backend():
41 have_any = False
42 for cmd in [[exe] for exe in program_bins.values()]:
43 try:
44 p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE)
45 (stdout, stderr) = p.communicate()
46 have_any = True
48 except OSError:
49 pass
51 return have_any
54qseis_components = 'r t z v'.split()
55qseis_greenf_names = ('ex', 'ss', 'ds', 'cl', 'fz', 'fh')
58def nextpow2(i):
59 return 2**int(math.ceil(math.log(i)/math.log(2.)))
62def str_float_vals(vals):
63 return ' '.join('%e' % val for val in vals)
66def str_int_vals(vals):
67 return ' '.join('%i' % val for val in vals)
70def str_str_vals(vals):
71 return ' '.join("'%s'" % val for val in vals)
74def scl(cs):
75 if not cs:
76 return '\n#'
78 return '\n'+' '.join('(%e,%e)' % (c.real, c.imag) for c in cs)
81def cake_model_to_config(mod):
82 k = 1000.
83 srows = []
84 ref_depth = 0.
85 for i, row in enumerate(mod.to_scanlines()):
86 depth, vp, vs, rho, qp, qs = row
87 if i == 0:
88 ref_depth = depth
90 row = [(depth-ref_depth)/k, vp/k, vs/k, rho/k, qp, qs]
91 srows.append('%i %s' % (i+1, str_float_vals(row)))
93 return '\n'.join(srows), len(srows), ref_depth
96class QSeisSourceMech(Object):
97 pass
100class QSeisSourceMechMT(QSeisSourceMech):
101 mnn = Float.T(default=1.0)
102 mee = Float.T(default=1.0)
103 mdd = Float.T(default=1.0)
104 mne = Float.T(default=0.0)
105 mnd = Float.T(default=0.0)
106 med = Float.T(default=0.0)
108 def string_for_config(self):
109 return '1 %(mnn)15e %(mee)15e %(mdd)15e ' \
110 '%(mne)15e %(med)15e %(mnd)15e ' % self.__dict__
113class QSeisSourceMechSDR(QSeisSourceMech):
114 m_iso = Float.T(default=0.0)
115 m_clvd = Float.T(default=0.0)
116 m_dc = Float.T(default=1.0e9)
117 strike = Float.T(default=0.0)
118 dip = Float.T(default=90.0)
119 rake = Float.T(default=0.0)
121 def string_for_config(self):
122 return '2 %(m_iso)15e %(m_clvd)15e %(m_dc)15e ' \
123 '%(strike)15e %(dip)15e %(rake)15e ' % self.__dict__
126class QSeisPropagationFilter(Object):
127 min_depth = Float.T(default=0.0)
128 max_depth = Float.T(default=0.0)
129 filtered_phase = Int.T(default=0)
131 def string_for_config(self):
132 return '%(min_depth)15e %(max_depth)15e ' \
133 '%(filtered_phase)i' % self.__dict__
136class QSeisPoleZeroFilter(Object):
137 constant = Complex.T(default=(1+0j))
138 poles = List.T(Complex.T())
139 zeros = List.T(Complex.T())
141 def string_for_config(self, version=None):
142 if version in ('2006a', '2006b'):
143 return '(%e,%e)\n%i%s\n%i%s' % (
144 self.constant.real, self.constant.imag,
145 len(self.zeros), scl(self.zeros),
146 len(self.poles), scl(self.poles))
147 elif version == '2006':
148 return '%e\n%i%s\n%i%s' % (
149 abs(self.constant),
150 len(self.zeros), scl(self.zeros),
151 len(self.poles), scl(self.poles))
154class QSeisConfig(Object):
156 qseis_version = String.T(default='2006')
157 time_region = Tuple.T(2, Timing.T(), default=(
158 Timing('-10'), Timing('+890')))
160 cut = Tuple.T(2, Timing.T(), optional=True)
161 fade = Tuple.T(4, Timing.T(), optional=True)
162 relevel_with_fade_in = Bool.T(default=False)
164 sw_algorithm = Int.T(default=0)
165 slowness_window = Tuple.T(4, Float.T(default=0.0))
166 wavenumber_sampling = Float.T(default=2.5)
167 aliasing_suppression_factor = Float.T(default=0.1)
168 source_disk_radius = Float.T(optional=True)
170 filter_surface_effects = Int.T(default=0)
171 filter_shallow_paths = Int.T(default=0)
172 filter_shallow_paths_depth = Float.T(default=0.0)
173 propagation_filters = List.T(QSeisPropagationFilter.T())
174 receiver_filter = QSeisPoleZeroFilter.T(optional=True)
176 sw_flat_earth_transform = Int.T(default=0)
178 gradient_resolution_vp = Float.T(default=0.0)
179 gradient_resolution_vs = Float.T(default=0.0)
180 gradient_resolution_density = Float.T(default=0.0)
182 wavelet_duration_samples = Float.T(default=0.0)
183 wavelet_type = Int.T(default=2)
184 user_wavelet_samples = List.T(Float.T())
186 def items(self):
187 return dict(self.T.inamevals(self))
190class QSeisConfigFull(QSeisConfig):
192 time_start = Float.T(default=0.0)
193 time_reduction_velocity = Float.T(default=0.0)
194 time_window = Float.T(default=900.0)
196 source_depth = Float.T(default=10.0)
197 source_mech = QSeisSourceMech.T(
198 optional=True,
199 default=QSeisSourceMechMT.D())
201 receiver_depth = Float.T(default=0.0)
202 receiver_distances = List.T(Float.T())
203 nsamples = Int.T(default=256)
205 gf_sw_source_types = Tuple.T(6, Int.T(), default=(1, 1, 1, 1, 0, 0))
207 gf_filenames = Tuple.T(6, String.T(), default=qseis_greenf_names)
209 seismogram_filename = String.T(default='seis')
211 receiver_azimuths = List.T(Float.T())
213 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True)
214 earthmodel_receiver_1d = gf.meta.Earthmodel1D.T(optional=True)
216 @staticmethod
217 def example():
218 conf = QSeisConfigFull()
219 conf.receiver_distances = [2000.]
220 conf.receiver_azimuths = [0.]
221 conf.time_start = -10.0
222 conf.time_reduction_velocity = 15.0
223 conf.earthmodel_1d = cake.load_model().extract(depth_max='cmb')
224 conf.earthmodel_receiver_1d = None
225 conf.sw_flat_earth_transform = 1
226 return conf
228 def get_output_filenames(self, rundir):
229 return [pjoin(rundir, self.seismogram_filename+'.t'+c)
230 for c in qseis_components]
232 def get_output_filenames_gf(self, rundir):
233 return [pjoin(rundir, fn+'.t'+c)
234 for fn in self.gf_filenames for c in qseis_components]
236 def string_for_config(self):
238 def aggregate(xx):
239 return len(xx), '\n'.join(
240 [''] + [x.string_for_config() for x in xx])
242 assert len(self.receiver_distances) > 0
243 assert len(self.receiver_distances) == len(self.receiver_azimuths)
244 assert self.earthmodel_1d is not None
246 d = self.__dict__.copy()
248 # fixing these switches here to reduce the amount of wrapper code
249 d['sw_distance_unit'] = 1 # always give distances in [km]
250 d['sw_t_reduce'] = 1 # time reduction always as velocity [km/s]
251 d['sw_equidistant'] = 0 # always give all distances and azimuths
252 d['sw_irregular_azimuths'] = 1
254 d['n_distances'] = len(self.receiver_distances)
255 d['str_distances'] = str_float_vals(self.receiver_distances)
256 d['str_azimuths'] = str_float_vals(self.receiver_azimuths)
258 model_str, nlines, ref_depth = cake_model_to_config(self.earthmodel_1d)
259 d['n_model_lines'] = nlines
260 d['model_lines'] = model_str
262 if self.earthmodel_receiver_1d:
263 model_str, nlines, ref_depth2 = cake_model_to_config(
264 self.earthmodel_receiver_1d)
265 assert ref_depth == ref_depth2
266 else:
267 model_str = "# no receiver side model"
268 nlines = 0
270 d['n_model_receiver_lines'] = nlines
271 d['model_receiver_lines'] = model_str
273 d['str_slowness_window'] = str_float_vals(self.slowness_window)
275 if self.qseis_version == '2006b':
276 sdr = self.source_disk_radius
277 d['str_source_disk_radius'] \
278 = '\n %e |dble: source_radius;' % (
279 sdr if sdr is not None else -1.0)
280 else:
281 if self.source_disk_radius is not None:
282 raise QSeisError(
283 'This version of QSEIS does not support the '
284 '`source_disk_radius` parameter.')
286 d['str_source_disk_radius'] = ''
288 if self.propagation_filters and ref_depth != 0.0:
289 raise QSeisError(
290 'Earth model must start with zero depth if '
291 'propagation_filters are set.')
293 d['n_depth_ranges'], d['str_depth_ranges'] = \
294 aggregate(self.propagation_filters)
296 if self.wavelet_type == 0: # user wavelet
297 d['str_w_samples'] = '\n' \
298 + '%i\n' % len(self.user_wavelet_samples) \
299 + str_float_vals(self.user_wavelet_samples)
300 else:
301 d['str_w_samples'] = ''
303 if self.receiver_filter:
304 d['str_receiver_filter'] = self.receiver_filter.string_for_config(
305 self.qseis_version)
306 else:
307 if self.qseis_version in ('2006a', '2006b'):
308 d['str_receiver_filter'] = '(1.0,0.0)\n0\n#\n0'
309 else:
310 d['str_receiver_filter'] = '1.0\n0\n#\n0'
312 d['str_gf_sw_source_types'] = str_int_vals(self.gf_sw_source_types)
313 d['str_gf_filenames'] = str_str_vals(self.gf_filenames)
315 if self.source_mech:
316 d['str_source'] = '%s \'%s\'' % (
317 self.source_mech.string_for_config(),
318 self.seismogram_filename)
319 else:
320 d['str_source'] = '0'
322 d['source_depth_rel'] = d['source_depth'] - ref_depth / km
323 d['receiver_depth_rel'] = d['receiver_depth'] - ref_depth / km
325 template = '''# autogenerated QSEIS input by qseis.py
326#
327# This is the input file of FORTRAN77 program "qseis06" for calculation of
328# synthetic seismograms based on a layered halfspace earth model.
329#
330# by
331# Rongjiang Wang <wang@gfz-potsdam.de>
332# GeoForschungsZentrum Potsdam
333# Telegrafenberg, D-14473 Potsdam, Germany
334#
335# Last modified: Potsdam, Nov., 2006
336#
337# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
338# If not specified, SI Unit System is used overall!
339#
340# Coordinate systems:
341# cylindrical (z,r,t) with z = downward,
342# r = from source outward,
343# t = azmuth angle from north to east;
344# cartesian (x,y,z) with x = north,
345# y = east,
346# z = downward;
347# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
348#
349# SOURCE PARAMETERS
350# =================
351# 1. source depth [km]
352#------------------------------------------------------------------------------
353 %(source_depth_rel)e |dble: source_depth;
354#------------------------------------------------------------------------------
355#
356# RECEIVER PARAMETERS
357# ===================
358# 1. receiver depth [km]
359# 2. switch for distance sampling role (1/0 = equidistant/irregular); switch
360# for unit used (1/0 = km/deg)
361# 3. number of distance samples
362# 4. if equidistant, then start and end trace distance (> 0); else distance
363# list (please order the receiver distances from small to large)
364# 5. (reduced) time begin [sec] & length of time window [sec], number of time
365# samples (<= 2*nfmax in qsglobal.h)
366# 6. switch for unit of the following time reduction parameter: 1 = velocity
367# [km/sec], 0 = slowness [sec/deg]; time reduction parameter
368#------------------------------------------------------------------------------
369 %(receiver_depth_rel)e |dble: receiver_depth;
370 %(sw_equidistant)i %(sw_distance_unit)i |int: sw_equidistant, sw_d_unit;
371 %(n_distances)i |int: no_distances;
372 %(str_distances)s |dble: d_1,d_n; or d_1,d_2, ...(no comments in between!);
373 %(time_start)e %(time_window)e %(nsamples)i |dble: t_start,t_window; int: no_t_samples;
374 %(sw_t_reduce)i %(time_reduction_velocity)e |int: sw_t_reduce; dble: t_reduce;
375#------------------------------------------------------------------------------
376#
377# WAVENUMBER INTEGRATION PARAMETERS
378# =================================
379# 1. select slowness integration algorithm (0 = suggested for full wave-field
380# modelling; 1 or 2 = suggested when using a slowness window with narrow
381# taper range - a technique for suppressing space-domain aliasing);
382# 2. 4 parameters for low and high slowness (Note 1) cut-offs [s/km] with
383# tapering: 0 < slw1 < slw2 defining cosine taper at the lower end, and 0 <
384# slw3 < slw4 defining the cosine taper at the higher end. default values
385# will be used in case of inconsistent input of the cut-offs (possibly with
386# much more computational effort);
387# 3. parameter for sampling rate of the wavenumber integration (1 = sampled
388# with the spatial Nyquist frequency, 2 = sampled with twice higher than
389# the Nyquist, and so on: the larger this parameter, the smaller the space-
390# domain aliasing effect, but also the more computation effort);
391# 4. the factor for suppressing time domain aliasing (> 0 and <= 1) (Note 2).
392#------------------------------------------------------------------------------
393 %(sw_algorithm)i |int: sw_algorithm;
394 %(str_slowness_window)s |dble: slw(1-4);
395 %(wavenumber_sampling)e |dble: sample_rate;
396 %(aliasing_suppression_factor)e |dble: supp_factor;%(str_source_disk_radius)s
397#------------------------------------------------------------------------------
398#
399# OPTIONS FOR PARTIAL SOLUTIONS
400# (only applied to the source-site structure)
401# ===========================================
402#
403# 1. switch for filtering free surface effects (0 = with free surface, i.e.,
404# do not select this filter; 1 = without free surface; 2 = without free
405# surface but with correction on amplitude and wave form. Note switch 2
406# can only be used for receivers at the surface)
407# 2. switch for filtering waves with a shallow penetration depth (concerning
408# their whole trace from source to receiver), penetration depth limit [km]
409#
410# if this option is selected, waves whose travel path never exceeds the
411# given depth limit will be filtered ("seismic nuting"). the condition for
412# selecting this filter is that the given shallow path depth limit should
413# be larger than both source and receiver depth.
414#
415# 3. number of depth ranges where the following selected up/down-sp2oing P or
416# SV waves should be filtered
417# 4. the 1. depth range: upper and lower depth [km], switch for filtering P
418# or SV wave in this depth range:
419#
420# switch no: 1 2 3 4 other
421# filtered phase: P(up) P(down) SV(up) SV(down) Error
422#
423# 5. the 2. ...
424#
425# The partial solution options are useful tools to increase the numerical
426# significance of desired wave phases. Especially when the desired phases
427# are smaller than the undesired phases, these options should be selected
428# and carefully combined.
429#------------------------------------------------------------------------------
430 %(filter_surface_effects)i |int: isurf;
431 %(filter_shallow_paths)i %(filter_shallow_paths_depth)e |int: sw_path_filter; dble:shallow_depth_limit;
432 %(n_depth_ranges)i %(str_depth_ranges)s
433#------------------------------------------------------------------------------
434#
435# SOURCE TIME FUNCTION (WAVELET) PARAMETERS (Note 3)
436# ==================================================
437# 1. wavelet duration [unit = time sample rather than sec!], that is about
438# equal to the half-amplitude cut-off period of the wavelet (> 0. if <= 0,
439# then default value = 2 time samples will be used), and switch for the
440# wavelet form (0 = user's own wavelet; 1 = default wavelet: normalized
441# square half-sinusoid for simulating a physical delta impulse; 2 = tapered
442# Heaviside wavelet, i.e. integral of wavelet 1)
443# 2. IF user's own wavelet is selected, then number of the wavelet time samples
444# (<= 1024), and followed by
445# 3. equidistant wavelet time samples
446# 4 ...(continue) (! no comment lines allowed between the time sample list!)
447# IF default, delete line 2, 3, 4 ... or comment them out!
448#------------------------------------------------------------------------------
449 %(wavelet_duration_samples)e %(wavelet_type)i%(str_w_samples)s
450#------------------------------------------------------------------------------
451#
452# FILTER PARAMETERS OF RECEIVERS (SEISMOMETERS OR HYDROPHONES)
453# ============================================================
454# 1. constant coefficient (normalization factor)
455# 2. number of roots (<= nrootmax in qsglobal.h)
456# 3. list of the root positions in the complex format (Re,Im). If no roots,
457# comment out this line
458# 4. number of poles (<= npolemax in qsglobal.h)
459# 5. list of the pole positions in the complex format (Re,Im). If no poles,
460# comment out this line
461#------------------------------------------------------------------------------
462 %(str_receiver_filter)s
463#------------------------------------------------------------------------------
464#
465# OUTPUT FILES FOR GREEN'S FUNCTIONS (Note 4)
466# ===========================================
467# 1. selections of source types (yes/no = 1/0)
468# 2. file names of Green's functions (please give the names without extensions,
469# which will be appended by the program automatically: *.tz, *.tr, *.tt
470# and *.tv are for the vertical, radial, tangential, and volume change (for
471# hydrophones) components, respectively)
472#------------------------------------------------------------------------------
473# explosion strike-slip dip-slip clvd single_f_v single_f_h
474#------------------------------------------------------------------------------
475 %(str_gf_sw_source_types)s
476 %(str_gf_filenames)s
477#------------------------------------------------------------------------------
478# OUTPUT FILES FOR AN ARBITRARY POINT DISLOCATION SOURCE
479# (for applications to earthquakes)
480# ======================================================
481# 1. selection (0 = not selected; 1 or 2 = selected), if (selection = 1), then
482# the 6 moment tensor elements [N*m]: Mxx, Myy, Mzz, Mxy, Myz, Mzx (x is
483# northward, y is eastward and z is downard); else if (selection = 2), then
484# Mis [N*m] = isotropic moment part = (MT+MN+MP)/3, Mcl = CLVD moment part
485# = (2/3)(MT+MP-2*MN), Mdc = double-couple moment part = MT-MN, Strike [deg],
486# Dip [deg] and Rake [deg].
487#
488# Note: to use this option, the Green's functions above should be computed
489# (selection = 1) if they do not exist already.
490#
491# north(x)
492# /
493# /\ strike
494# *-----------------------> east(y)
495# |\ \
496# |-\ \
497# | \ fault plane \
498# |90 \ \
499# |-dip\ \
500# | \ \
501# | \ \
502# downward(z) \-----------------------\\
503#
504# 2. switch for azimuth distribution of the stations (0 = uniform azimuth,
505# else = irregular azimuth angles)
506# 3. list of the azimuth angles [deg] for all stations given above (if the
507# uniform azimuth is selected, then only one azimuth angle is required)
508#
509#------------------------------------------------------------------------------
510# Mis Mcl Mdc Strike Dip Rake File
511#------------------------------------------------------------------------------
512# 2 0.00 1.00 6.0E+19 120.0 30.0 25.0 'seis'
513#------------------------------------------------------------------------------
514# Mxx Myy Mzz Mxy Myz Mzx File
515#------------------------------------------------------------------------------
516%(str_source)s
517%(sw_irregular_azimuths)i
518%(str_azimuths)s
519#------------------------------------------------------------------------------
520#
521# GLOBAL MODEL PARAMETERS (Note 5)
522# ================================
523# 1. switch for flat-earth-transform
524# 2. gradient resolution [%%] of vp, vs, and ro (density), if <= 0, then default
525# values (depending on wave length at cut-off frequency) will be used
526#------------------------------------------------------------------------------
527 %(sw_flat_earth_transform)i |int: sw_flat_earth_transform;
528 %(gradient_resolution_vp)e %(gradient_resolution_vs)e %(gradient_resolution_density)e |dble: vp_res, vs_res, ro_res;
529#------------------------------------------------------------------------------
530#
531# LAYERED EARTH MODEL
532# (SHALLOW SOURCE + UNIFORM DEEP SOURCE/RECEIVER STRUCTURE)
533# =========================================================
534# 1. number of data lines of the layered model (source site)
535#------------------------------------------------------------------------------
536 %(n_model_lines)i |int: no_model_lines;
537#------------------------------------------------------------------------------
538#
539# MULTILAYERED MODEL PARAMETERS (source site)
540# ===========================================
541# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs
542#------------------------------------------------------------------------------
543%(model_lines)s
544#------------------------------------------------------------------------------
545#
546# LAYERED EARTH MODEL
547# (ONLY THE SHALLOW RECEIVER STRUCTURE)
548# =====================================
549# 1. number of data lines of the layered model
550#
551# Note: if the number = 0, then the receiver site is the same as the
552# source site, else different receiver-site structure is considered.
553# please be sure that the lowest interface of the receiver-site
554# structure given given below can be found within the source-site
555# structure, too.
556#
557#------------------------------------------------------------------------------
558 %(n_model_receiver_lines)i |int: no_model_lines;
559#------------------------------------------------------------------------------
560#
561# MULTILAYERED MODEL PARAMETERS (shallow receiver-site structure)
562# ===============================================================
563# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs
564#------------------------------------------------------------------------------
565%(model_receiver_lines)s
566#---------------------------------end of all inputs----------------------------
569Note 1:
571The slowness is defined by inverse value of apparent wave velocity = sin(i)/v
572with i = incident angle and v = true wave velocity.
574Note 2:
576The suppression of the time domain aliasing is achieved by using the complex
577frequency technique. The suppression factor should be a value between 0 and 1.
578If this factor is set to 0.1, for example, the aliasing phase at the reduced
579time begin is suppressed to 10%%.
581Note 3:
583The default basic wavelet function (option 1) is (2/tau)*sin^2(pi*t/tau),
584for 0 < t < tau, simulating physical delta impuls. Its half-amplitude cut-off
585frequency is 1/tau. To avoid high-frequency noise, tau should not be smaller
586than 4-5 time samples.
588Note 4:
590 Double-Couple m11/ m22/ m33/ m12/ m23/ m31 Azimuth_Factor_(tz,tr,tv)/(tt)
591 ============================================================================
592 explosion 1.0/ 1.0/ 1.0/ -- / -- / -- 1.0 / 0.0
593 strike-slip -- / -- / -- / 1.0/ -- / -- sin(2*azi) / cos(2*azi)
594 1.0/-1.0/ -- / -- / -- / -- cos(2*azi) / -sin(2*azi)
595 dip-slip -- / -- / -- / -- / -- / 1.0 cos(azi) / sin(azi)
596 -- / -- / -- / -- / 1.0/ -- sin(azi) / -cos(azi)
597 clvd -0.5/-0.5/ 1.0/ -- / -- / -- 1.0 / 0.0
598 ============================================================================
599 Single-Force fx / fy / fz Azimuth_Factor_(tz,tr,tv)/(tt)
600 ============================================================================
601 fz -- / -- / 1.0 1.0 / 0.0
602 fx 1.0/ -- / -- cos(azi) / sin(azi)
603 fy -- / 1.0/ -- sin(azi) / -cos(azi)
604 ============================================================================
606Note 5:
608Layers with a constant gradient will be discretized with a number of homogeneous
609sublayers. The gradient resolutions are then used to determine the maximum
610allowed thickness of the sublayers. If the resolutions of Vp, Vs and Rho
611(density) require different thicknesses, the smallest is first chosen. If this
612is even smaller than 1%% of the characteristic wavelength, then the latter is
613taken finally for the sublayer thickness.
614''' # noqa
616 return (template % d).encode('ascii')
619class QSeisError(gf.store.StoreError):
620 pass
623class Interrupted(gf.store.StoreError):
624 def __str__(self):
625 return 'Interrupted.'
628class QSeisRunner(object):
630 def __init__(self, tmp=None, keep_tmp=False):
631 self.tempdir = mkdtemp(prefix='qseisrun-', dir=tmp)
632 self.keep_tmp = keep_tmp
633 self.config = None
635 def run(self, config):
636 self.config = config
638 input_fn = pjoin(self.tempdir, 'input')
640 with open(input_fn, 'wb') as f:
641 input_str = config.string_for_config()
642 logger.debug('===== begin qseis input =====\n'
643 '%s===== end qseis input =====' % input_str.decode())
644 f.write(input_str)
646 program = program_bins['qseis.%s' % config.qseis_version]
648 old_wd = os.getcwd()
650 os.chdir(self.tempdir)
652 interrupted = []
654 def signal_handler(signum, frame):
655 os.kill(proc.pid, signal.SIGTERM)
656 interrupted.append(True)
658 original = signal.signal(signal.SIGINT, signal_handler)
659 try:
660 try:
661 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
663 except OSError:
664 os.chdir(old_wd)
665 raise QSeisError(
666 '''could not start qseis executable: "%s"
667Available fomosto backends and download links to the modelling codes are listed
668on
670 https://pyrocko.org/docs/current/apps/fomosto/backends.html
672''' % program)
674 (output_str, error_str) = proc.communicate(b'input\n')
676 finally:
677 signal.signal(signal.SIGINT, original)
679 if interrupted:
680 raise KeyboardInterrupt()
682 logger.debug('===== begin qseis output =====\n'
683 '%s===== end qseis output =====' % output_str.decode())
685 errmess = []
686 if proc.returncode != 0:
687 errmess.append(
688 'qseis had a non-zero exit state: %i' % proc.returncode)
690 if error_str:
691 logger.warning(
692 'qseis emitted something via stderr:\n\n%s'
693 % error_str.decode())
695 # errmess.append('qseis emitted something via stderr')
697 if output_str.lower().find(b'error') != -1:
698 errmess.append("the string 'error' appeared in qseis output")
700 if errmess:
701 self.keep_tmp = True
703 os.chdir(old_wd)
704 raise QSeisError('''
705===== begin qseis input =====
706%s===== end qseis input =====
707===== begin qseis output =====
708%s===== end qseis output =====
709===== begin qseis error =====
710%s===== end qseis error =====
711%s
712qseis has been invoked as "%s"
713in the directory %s'''.lstrip() % (
714 input_str.decode(), output_str.decode(), error_str.decode(),
715 '\n'.join(errmess), program, self.tempdir))
717 self.qseis_output = output_str
718 self.qseis_error = error_str
720 os.chdir(old_wd)
722 def get_traces(self, which='seis'):
724 if which == 'seis':
725 fns = self.config.get_output_filenames(self.tempdir)
726 components = qseis_components
728 elif which == 'gf':
729 fns = self.config.get_output_filenames_gf(self.tempdir)
730 components = [
731 fn+'.t'+c
732 for fn in self.config.gf_filenames for c in qseis_components]
733 else:
734 raise Exception(
735 'get_traces: which argument should be "seis" or "gf"')
737 traces = []
738 distances = self.config.receiver_distances
739 azimuths = self.config.receiver_azimuths
740 for comp, fn in zip(components, fns):
741 if not os.path.exists(fn):
742 continue
744 data = num.loadtxt(fn, skiprows=1, dtype=float)
745 nsamples, ntraces = data.shape
746 ntraces -= 1
747 vred = self.config.time_reduction_velocity
748 deltat = (data[-1, 0] - data[0, 0])/(nsamples-1)
750 for itrace, distance, azimuth in zip(
751 range(ntraces), distances, azimuths):
753 tmin = self.config.time_start
754 if vred != 0.0:
755 tmin += distance / vred
757 tmin += deltat
758 tr = trace.Trace(
759 '', '%04i' % itrace, '', comp,
760 tmin=tmin, deltat=deltat, ydata=data[:, itrace+1],
761 meta=dict(
762 distance=distance*km,
763 azimuth=azimuth))
765 traces.append(tr)
767 return traces
769 def __del__(self):
770 if self.tempdir:
771 if not self.keep_tmp:
772 shutil.rmtree(self.tempdir)
773 self.tempdir = None
774 else:
775 logger.warning(
776 'not removing temporary directory: %s' % self.tempdir)
779class QSeisGFBuilder(gf.builder.Builder):
780 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
781 force=False):
783 self.store = gf.store.Store(store_dir, 'w')
785 if block_size is None:
786 block_size = (1, 1, 100)
788 if len(self.store.config.ns) == 2:
789 block_size = block_size[1:]
791 gf.builder.Builder.__init__(
792 self, self.store.config, step, block_size=block_size, force=force)
794 baseconf = self.store.get_extra('qseis')
796 conf = QSeisConfigFull(**baseconf.items())
797 conf.earthmodel_1d = self.store.config.earthmodel_1d
798 conf.earthmodel_receiver_1d = self.store.config.earthmodel_receiver_1d
800 deltat = 1.0/self.gf_config.sample_rate
802 if 'time_window_min' not in shared:
803 d = self.store.make_timing_params(
804 conf.time_region[0], conf.time_region[1],
805 force=force)
807 shared['time_window_min'] = d['tlenmax_vred']
808 shared['time_start'] = d['tmin_vred']
809 shared['time_reduction_velocity'] = d['vred'] / km
811 time_window_min = shared['time_window_min']
812 conf.time_start = shared['time_start']
814 conf.time_reduction_velocity = shared['time_reduction_velocity']
816 conf.nsamples = nextpow2(int(round(time_window_min / deltat)) + 1)
817 conf.time_window = (conf.nsamples-1)*deltat
819 self.qseis_config = conf
821 self.tmp = tmp
822 if self.tmp is not None:
823 util.ensuredir(self.tmp)
825 def cleanup(self):
826 self.store.close()
828 def work_block(self, index):
829 if len(self.store.config.ns) == 2:
830 (sz, firstx), (sz, lastx), (ns, nx) = \
831 self.get_block_extents(index)
833 rz = self.store.config.receiver_depth
834 else:
835 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
836 self.get_block_extents(index)
838 conf = copy.deepcopy(self.qseis_config)
840 logger.info('Starting block %i / %i' %
841 (index+1, self.nblocks))
843 conf.source_depth = float(sz/km)
844 conf.receiver_depth = float(rz/km)
846 runner = QSeisRunner(tmp=self.tmp)
848 dx = self.gf_config.distance_delta
850 distances = num.linspace(firstx, firstx + (nx-1)*dx, nx).tolist()
852 if distances[-1] < self.gf_config.distance_max:
853 # add global max distance, because qseis does some adjustments with
854 # this value
855 distances.append(self.gf_config.distance_max)
857 mex = (MomentTensor(m=symmat6(1, 1, 1, 0, 0, 0)),
858 {'r': (0, +1), 'z': (1, +1)})
860 mmt1 = (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
861 {'r': (0, +1), 't': (3, +1), 'z': (5, +1)})
862 mmt2 = (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
863 {'r': (1, +1), 't': (4, +1), 'z': (6, +1)})
864 mmt3 = (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
865 {'r': (2, +1), 'z': (7, +1)})
866 mmt4 = (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
867 {'r': (8, +1), 'z': (9, +1)})
869 component_scheme = self.store.config.component_scheme
870 off = 0
871 if component_scheme == 'elastic8':
872 off = 8
873 elif component_scheme == 'elastic10':
874 off = 10
876 msf = (None, {
877 'fz.tr': (off+0, +1),
878 'fh.tr': (off+1, +1),
879 'fh.tt': (off+2, -1),
880 'fz.tz': (off+3, +1),
881 'fh.tz': (off+4, +1)})
882 if component_scheme == 'elastic2':
883 gfsneeded = (1, 0, 0, 0, 0, 0)
884 gfmapping = [mex]
886 if component_scheme == 'elastic5':
887 gfsneeded = (0, 0, 0, 0, 1, 1)
888 gfmapping = [msf]
890 elif component_scheme == 'elastic8':
891 gfsneeded = (1, 1, 1, 1, 0, 0)
892 gfmapping = [mmt1, mmt2, mmt3]
894 elif component_scheme == 'elastic10':
895 gfsneeded = (1, 1, 1, 1, 0, 0)
896 gfmapping = [mmt1, mmt2, mmt3, mmt4]
898 elif component_scheme == 'elastic13':
899 gfsneeded = (1, 1, 1, 1, 1, 1)
900 gfmapping = [mmt1, mmt2, mmt3, msf]
902 elif component_scheme == 'elastic15':
903 gfsneeded = (1, 1, 1, 1, 1, 1)
904 gfmapping = [mmt1, mmt2, mmt3, mmt4, msf]
906 conf.gf_sw_source_types = gfsneeded
907 conf.receiver_distances = [d/km for d in distances]
908 conf.receiver_azimuths = [0.0] * len(distances)
910 for mt, gfmap in gfmapping:
911 if mt:
912 m = mt.m()
913 f = float
914 conf.source_mech = QSeisSourceMechMT(
915 mnn=f(m[0, 0]), mee=f(m[1, 1]), mdd=f(m[2, 2]),
916 mne=f(m[0, 1]), mnd=f(m[0, 2]), med=f(m[1, 2]))
917 else:
918 conf.source_mech = None
920 if any(conf.gf_sw_source_types) or conf.source_mech is not None:
921 runner.run(conf)
923 if any(c in gfmap for c in qseis_components):
924 rawtraces = runner.get_traces('seis')
925 else:
926 rawtraces = runner.get_traces('gf')
928 interrupted = []
930 def signal_handler(signum, frame):
931 interrupted.append(True)
933 original = signal.signal(signal.SIGINT, signal_handler)
934 self.store.lock()
935 try:
936 for itr, tr in enumerate(rawtraces):
937 if tr.channel not in gfmap:
938 continue
940 x = tr.meta['distance']
941 if x > firstx + (nx-1)*dx:
942 continue
944 ig, factor = gfmap[tr.channel]
946 if len(self.store.config.ns) == 2:
947 args = (sz, x, ig)
948 else:
949 args = (rz, sz, x, ig)
951 if conf.cut:
952 tmin = self.store.t(conf.cut[0], args[:-1])
953 tmax = self.store.t(conf.cut[1], args[:-1])
955 if None in (tmin, tmax):
956 self.warn(
957 'Failed cutting {} traces. ' +
958 'Failed to determine time window')
959 continue
961 tr.chop(tmin, tmax)
963 tmin = tr.tmin
964 tmax = tr.tmax
966 if conf.fade:
967 ta, tb, tc, td = [
968 self.store.t(v, args[:-1]) for v in conf.fade]
970 if None in (ta, tb, tc, td):
971 continue
973 if not (ta <= tb and tb <= tc and tc <= td):
974 raise QSeisError(
975 'invalid fade configuration '
976 '(it should be (ta <= tb <= tc <= td) but '
977 'ta=%g, tb=%g, tc=%g, td=%g)' % (
978 ta, tb, tc, td))
980 t = tr.get_xdata()
981 fin = num.interp(t, [ta, tb], [0., 1.])
982 fout = num.interp(t, [tc, td], [1., 0.])
983 anti_fin = 1. - fin
984 anti_fout = 1. - fout
986 y = tr.ydata
988 sum_anti_fin = num.sum(anti_fin)
989 sum_anti_fout = num.sum(anti_fout)
991 if sum_anti_fin != 0.0:
992 yin = num.sum(anti_fin*y) / sum_anti_fin
993 else:
994 yin = 0.0
996 if sum_anti_fout != 0.0:
997 yout = num.sum(anti_fout*y) / sum_anti_fout
998 else:
999 yout = 0.0
1001 y2 = anti_fin*yin + fin*fout*y + anti_fout*yout
1003 if conf.relevel_with_fade_in:
1004 y2 -= yin
1006 tr.set_ydata(y2)
1008 gf_tr = gf.store.GFTrace.from_trace(tr)
1009 gf_tr.data *= factor
1011 try:
1012 self.store.put(args, gf_tr)
1013 except gf.store.DuplicateInsert:
1014 self.warn('{} insertions_skipped (duplicates)')
1016 finally:
1017 self.log_warnings(index+1, logger)
1018 self.store.unlock()
1019 signal.signal(signal.SIGINT, original)
1021 if interrupted:
1022 raise KeyboardInterrupt()
1024 conf.gf_sw_source_types = (0, 0, 0, 0, 0, 0)
1026 logger.info('Done with block %i / %i' %
1027 (index+1, self.nblocks))
1030def init(store_dir, variant, config_params=None):
1031 if variant is None:
1032 variant = '2006'
1034 if variant not in ('2006', '2006a', '2006b'):
1035 raise gf.store.StoreError('unsupported variant: %s' % variant)
1037 modelling_code_id = 'qseis.%s' % variant
1039 qseis = QSeisConfig(qseis_version=variant)
1040 qseis.time_region = (
1041 gf.meta.Timing('begin-50'),
1042 gf.meta.Timing('end+100'))
1044 qseis.cut = (
1045 gf.meta.Timing('begin-50'),
1046 gf.meta.Timing('end+100'))
1048 qseis.wavelet_duration_samples = 0.001
1049 qseis.sw_flat_earth_transform = 1
1051 store_id = os.path.basename(os.path.realpath(store_dir))
1053 d = dict(
1054 id=store_id,
1055 ncomponents=10,
1056 sample_rate=0.2,
1057 receiver_depth=0*km,
1058 source_depth_min=10*km,
1059 source_depth_max=20*km,
1060 source_depth_delta=10*km,
1061 distance_min=100*km,
1062 distance_max=1000*km,
1063 distance_delta=10*km,
1064 earthmodel_1d=cake.load_model().extract(depth_max='cmb'),
1065 modelling_code_id=modelling_code_id,
1066 tabulated_phases=[
1067 gf.meta.TPDef(
1068 id='begin',
1069 definition='p,P,p\\,P\\,Pv_(cmb)p'),
1070 gf.meta.TPDef(
1071 id='end',
1072 definition='2.5'),
1073 gf.meta.TPDef(
1074 id='P',
1075 definition='!P'),
1076 gf.meta.TPDef(
1077 id='S',
1078 definition='!S'),
1079 gf.meta.TPDef(
1080 id='p',
1081 definition='!p'),
1082 gf.meta.TPDef(
1083 id='s',
1084 definition='!s')])
1086 if config_params is not None:
1087 d.update(config_params)
1089 config = gf.meta.ConfigTypeA(**d)
1091 config.validate()
1092 return gf.store.Store.create_editables(
1093 store_dir, config=config, extra={'qseis': qseis})
1096def build(store_dir, force=False, nworkers=None, continue_=False, step=None,
1097 iblock=None):
1099 return QSeisGFBuilder.build(
1100 store_dir, force=force, nworkers=nworkers, continue_=continue_,
1101 step=step, iblock=iblock)