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