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