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