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
17import os.path as op
18from scipy.integrate import cumtrapz
20from pyrocko.moment_tensor import MomentTensor, symmat6
21from pyrocko.guts import Float, Int, Tuple, List, Bool, Object, String
22from pyrocko import trace, util, cake, gf
24km = 1e3
26guts_prefix = 'pf'
28Timing = gf.meta.Timing
31logger = logging.getLogger('pyrocko.fomosto.qseis2d')
33# how to call the programs
34program_bins = {
35 'qseis2d.qseisS2014': 'fomosto_qseisS2014',
36 'qseis2d.qseisR2014': 'fomosto_qseisR2014',
37}
39qseis2d_components = {
40 1: 'z r t'.split(),
41 2: 'e n u'.split(),
42}
44# defaults
45default_gf_directory = 'qseis2d_green'
46default_fk_basefilename = 'green'
47default_source_depth = 10.0
48default_time_region = (Timing('-10'), Timing('+890'))
49default_slowness_window = (0.0, 0.0, 0.2, 0.25)
52def have_backend():
53 for cmd in [[exe] for exe in program_bins.values()]:
54 try:
55 p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE)
56 (stdout, stderr) = p.communicate()
58 except OSError:
59 return False
61 return True
64def nextpow2(i):
65 return 2 ** int(math.ceil(math.log(i) / math.log(2.)))
68def str_float_vals(vals):
69 return ' '.join('%e' % val for val in vals)
72def str_int_vals(vals):
73 return ' '.join('%i' % val for val in vals)
76def str_str_vals(vals):
77 return ' '.join("'%s'" % val for val in vals)
80def scl(cs):
81 if not cs:
82 return '\n#'
84 return '\n' + ' '.join('(%e,%e)' % (c.real, c.imag) for c in cs)
87def cake_model_to_config(mod):
88 k = 1000.
89 srows = []
90 for i, row in enumerate(mod.to_scanlines()):
91 depth, vp, vs, rho, qp, qs = row
92 row = [depth / k, vp / k, vs / k, rho / k, qp, qs]
93 srows.append('%i %s' % (i + 1, str_float_vals(row)))
95 return '\n'.join(srows), len(srows)
98class QSeis2dSource(Object):
99 lat = Float.T(default=10.)
100 lon = Float.T(default=0.)
101 depth = Float.T(default=default_source_depth)
103 def string_for_config(self):
104 return '%(lat)e %(lon)15e ' % self.__dict__
107class QSeisRSourceMech(Object):
108 pass
111class QSeisRSourceMechMT(QSeisRSourceMech):
112 mnn = Float.T(default=1.0)
113 mee = Float.T(default=1.0)
114 mdd = Float.T(default=1.0)
115 mne = Float.T(default=0.0)
116 mnd = Float.T(default=0.0)
117 med = Float.T(default=0.0)
119 def string_for_config(self):
120 return '%(mnn)e %(mee)15e %(mdd)15e ' \
121 '%(mne)15e %(med)15e %(mnd)15e ' % self.__dict__
124class QSeisSPropagationFilter(Object):
125 min_depth = Float.T(default=0.0)
126 max_depth = Float.T(default=0.0)
127 filtered_phase = Int.T(default=0)
129 def string_for_config(self):
130 return '%(min_depth)15e %(max_depth)15e ' \
131 '%(filtered_phase)i' % self.__dict__
134class QSeisRBandpassFilter(Object):
135 order = Int.T(default=-1)
136 lower_cutoff = Float.T(default=0.1) # [Hz]
137 upper_cutoff = Float.T(default=10.0)
139 def string_for_config(self):
140 return '%(order) %(lower_cutoff)5e %(upper_cutoff)5e' % self.__dict__
143class QSeisSConfig(Object):
145 qseiss_version = String.T(default='2014')
147 calc_slowness_window = Int.T(default=1)
148 slowness_window = Tuple.T(4, optional=True)
149 wavenumber_sampling = Float.T(default=2.5)
150 aliasing_suppression_factor = Float.T(default=0.01)
152 filter_shallow_paths = Int.T(default=0)
153 filter_shallow_paths_depth = Float.T(default=0.0)
154 propagation_filters = List.T(QSeisSPropagationFilter.T())
156 sw_flat_earth_transform = Int.T(default=0)
158 gradient_resolution_vp = Float.T(default=0.0)
159 gradient_resolution_vs = Float.T(default=0.0)
160 gradient_resolution_density = Float.T(default=0.0)
162 def items(self):
163 return dict(self.T.inamevals(self))
166class QSeisSConfigFull(QSeisSConfig):
168 time_window = Float.T(default=900.0)
170 source_depth = Float.T(default=10.0)
172 receiver_basement_depth = Float.T(default=35.0) # [km]
173 receiver_min_distance = Float.T(default=1000.0) # [km]
174 receiver_max_distance = Float.T(default=10000.0) # [km]
175 nsamples = Int.T(default=256)
177 info_path = String.T(default='green.info')
178 fk_path = String.T(default='green.fk')
180 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True)
182 @staticmethod
183 def example():
184 conf = QSeisSConfigFull()
185 conf.source_depth = 15.
186 conf.receiver_basement_depth = 35.
187 conf.receiver_max_distance = 2000.
188 conf.earthmodel_1d = cake.load_model().extract(depth_max='cmb')
189 conf.sw_flat_earth_transform = 1
190 return conf
192 def string_for_config(self):
194 def aggregate(xx):
195 return len(xx), '\n'.join(
196 [''] + [x.string_for_config() for x in xx])
198 assert self.earthmodel_1d is not None
199 assert self.slowness_window is not None or self.calc_slowness_window,\
200 "'slowness window' undefined and 'calc_slowness_window'=0"
202 d = self.__dict__.copy()
204 model_str, nlines = cake_model_to_config(self.earthmodel_1d)
205 d['n_model_lines'] = nlines
206 d['model_lines'] = model_str
208 if not self.slowness_window:
209 d['str_slowness_window'] = str_float_vals(default_slowness_window)
210 else:
211 d['str_slowness_window'] = str_float_vals(self.slowness_window)
213 d['n_depth_ranges'], d['str_depth_ranges'] = \
214 aggregate(self.propagation_filters)
216 template = '''# autogenerated QSEISS input by qseis2d.py
217# This is the input file of FORTRAN77 program "QseisS" for calculation of
218# f-k spectra of upgoing seismic waves at the reveiver-site basement.
219#
220# by
221# Rongjiang Wang <wang@gfz-potsdam.de>
222# GeoForschungsZentrum Potsdam
223# Telegrafenberg, D-14473 Potsdam, Germany
224#
225# Modified from qseis2006, Potsdam, Dec. 2014
226#
227# SOURCE DEPTH
228# ============
229# 1. source depth [km]
230#------------------------------------------------------------------------------
231 %(source_depth)e |dble;
232#------------------------------------------------------------------------------
233#
234# RECEIVER-SITE PARAMETERS
235# ========================
236# 1. receiver-site basement depth [km]
237# 2. max. epicental distance [km]
238#------------------------------------------------------------------------------
239 %(receiver_basement_depth)e |dble;
240 %(receiver_max_distance)e |dble;
241#------------------------------------------------------------------------------
242# TIME SAMPLING PARAMETERS
243# ========================
244# 1. length of time window [sec]
245# 2. number of time samples (<= 2*nfmax in qsglobal.h)
246#------------------------------------------------------------------------------
247 %(time_window)e |dble: t_window;
248 %(nsamples)i |int: no_t_samples;
249#------------------------------------------------------------------------------
250# SLOWNESS WINDOW PARAMETERS
251# ==========================
252# 1. the low and high slowness cut-offs [s/km] with tapering: 0 < slw1 < slw2
253# defining the bounds of cosine taper at the lower end, and 0 < slw3 < slw4
254# defining the bounds of cosine taper at the higher end.
255# 2. slowness sampling rate (1.0 = sampling with the Nyquist slowness, 2.0 =
256# sampling with twice higher than the Nyquist slowness, and so on: the
257# larger this parameter, the smaller the space-domain aliasing effect, but
258# also the more computation effort);
259# 3. the factor for suppressing time domain aliasing (> 0 and <= 1).
260#------------------------------------------------------------------------------
261 %(str_slowness_window)s |dble: slw(1-4);
262 %(wavenumber_sampling)e |dble: sample_rate;
263 %(aliasing_suppression_factor)e |dble: supp_factor;
264#------------------------------------------------------------------------------
265# OPTIONS FOR PARTIAL SOLUTIONS
266# =============================
267# 1. switch for filtering waves with a shallow penetration depth (concerning
268# their whole trace from source to receiver), penetration depth limit [km]
269# (Note: if this option is selected, waves whose travel path never exceeds
270# the given depth limit will be filtered ("seismic nuting"). the condition
271# for selecting this filter is that the given shallow path depth limit
272# should be larger than both source and receiver depth.)
273# 2. number of depth ranges where the following selected up/down-sp2oing P or
274# SV waves should be filtered
275# 3. the 1. depth range: upper and lower depth [km], switch for filtering P
276# or SV wave in this depth range:
277# switch no: 1 2 3 4 other
278# filtered phase: P(up) P(down) SV(up) SV(down) Error
279# 4. the 2. ...
280# (Note: the partial solution options are useful tools to increase the
281# numerical resolution of desired wave phases. especially when the desired
282# phases are much smaller than the undesired phases, these options should
283# be selected and carefully combined.)
284#------------------------------------------------------------------------------
285 %(filter_shallow_paths)i %(filter_shallow_paths_depth)e
286 %(n_depth_ranges)i %(str_depth_ranges)s |int, dble, dble, int;
287#------------------------------------------------------------------------------
288# OUTPUT FILE FOR F-K SPECTRA of GREEN'S FUNCTIONS
289# ================================================
290# 1. info-file of Green's functions (ascii including f-k sampling parameters)
291# 2. file name of Green's functions (binary files including explosion, strike
292# -slip, dip-slip and clvd sources)
293#------------------------------------------------------------------------------
294 '%(info_path)s'
295 '%(fk_path)s'
296#------------------------------------------------------------------------------
297# GLOBAL MODEL PARAMETERS
298# =======================
299# 1. switch for flat-earth-transform
300# 2. gradient resolution [percent] of vp, vs, and rho (density), if <= 0, then
301# default values (depending on wave length at cut-off frequency) will be
302# used
303#------------------------------------------------------------------------------
304 %(sw_flat_earth_transform)i |int: sw_flat_earth_transform;
305 %(gradient_resolution_vp)e %(gradient_resolution_vs)e %(gradient_resolution_density)e
306#------------------------------------------------------------------------------
307# SOURCE-SITE LAYERED EARTH MODEL
308# ===============================
309# 1. number of data lines of the layered model
310#------------------------------------------------------------------------------
311 %(n_model_lines)i |int: no_model_lines;
312#------------------------------------------------------------------------------
313# no depth[km] vp[km/s] vs[km/s] rho[g/cm^3] qp qs
314#------------------------------------------------------------------------------
315%(model_lines)s ;
316#-----------------------END OF INPUT PARAMETERS--------------------------------
318Glossary
320SLOWNESS: The slowness is the inverse of apparent wave velocity = sin(i)/v with
321i = incident angle and v = true wave velocity.
323SUPPRESSION OF TIME-DOMAIN ALIASING: The suppression of the time domain aliasing
324is achieved by using the complex frequency technique. The suppression factor
325should be a value between 0 and 1. If this factor is set to 0.1, for example, the
326aliasing phase at the reduced time begin is suppressed to 10 percent.
328MODEL PARAMETER GRADIENT RESOLUTION: Layers with a constant gradient will be
329discretized with a number of homogeneous sublayers. The gradient resolutions are
330then used to determine the maximum allowed thickness of the sublayers. If the
331resolutions of Vp, Vs and Rho (density) require different thicknesses, the
332smallest is first chosen. If this is even smaller than 1 percent of the
333characteristic wavelength, then the latter is taken finally for the sublayer
334thickness.
335''' # noqa
336 return (template % d).encode('ascii')
339class QSeisRReceiver(Object):
340 lat = Float.T(default=10.0)
341 lon = Float.T(default=0.0)
342 depth = Float.T(default=0.0)
343 tstart = Float.T(default=0.0)
344 distance = Float.T(default=0.0)
346 def string_for_config(self):
347 return '%(lat)e %(lon)15e %(depth)15e' % self.__dict__
350class QSeisRConfig(Object):
352 qseisr_version = String.T(default='2014')
354 receiver_filter = QSeisRBandpassFilter.T(optional=True)
356 wavelet_duration = Float.T(default=0.001) # [s]
357 wavelet_type = Int.T(default=1)
358 user_wavelet_samples = List.T(Float.T())
360 def items(self):
361 return dict(self.T.inamevals(self))
364class QSeisRConfigFull(QSeisRConfig):
366 source = QSeis2dSource(depth=default_source_depth) # [lat, lon(deg)]
367 receiver = QSeisRReceiver() # [lat, lon(deg)]
369 source_mech = QSeisRSourceMech.T(
370 optional=True,
371 default=QSeisRSourceMechMT.D())
373 time_reduction = Float.T(default=0.0)
375 info_path = String.T(default='green.info')
376 fk_path = String.T(default='green.fk')
378 output_format = Int.T(default=1) # 1/2 components in [Z, R, T]/[E, N, U]
379 output_filename = String.T(default='seis.dat')
381 earthmodel_receiver_1d = gf.meta.Earthmodel1D.T(optional=True)
383 @staticmethod
384 def example():
385 conf = QSeisRConfigFull()
386 conf.source = QSeis2dSource(lat=-80.5, lon=90.1)
387 conf.receiver_location = QSeisRReceiver(lat=13.4, lon=240.5, depth=0.0)
388 conf.time_reduction = 10.0
389 conf.earthmodel_receiver_1d = cake.load_model().extract(
390 depth_max='moho')
391 return conf
393 @property
394 def components(self):
395 return qseis2d_components[self.output_format]
397 def get_output_filename(self, rundir):
398 return op.join(rundir, self.output_filename)
400 def string_for_config(self):
402 def aggregate(xx):
403 return len(xx), '\n'.join(
404 [''] + [x.string_for_config() for x in xx])
406 assert self.earthmodel_receiver_1d is not None
408 d = self.__dict__.copy()
410# Actually not existing anymore in code
411# d['sw_surface'] = 0 # 0-free-surface, 1-no fre surface
413 d['str_source_location'] = self.source.string_for_config()
415 d['str_receiver'] = self.receiver.string_for_config()
417 d['str_output_filename'] = "'%s'" % self.output_filename
419 model_str, nlines = cake_model_to_config(self.earthmodel_receiver_1d)
421 d['n_model_receiver_lines'] = nlines
422 d['model_receiver_lines'] = model_str
424 if self.wavelet_type == 0: # user wavelet
425 d['str_w_samples'] = '\n' \
426 + '%i\n' % len(self.user_wavelet_samples) \
427 + str_float_vals(self.user_wavelet_samples)
428 else:
429 d['str_w_samples'] = ''
431 if self.receiver_filter:
432 d['str_receiver_filter'] = self.receiver_filter.string_for_config()
433 else:
434 d['str_receiver_filter'] = '-1 0.0 200.'
436 if self.source_mech:
437 d['str_source'] = '%s' % (self.source_mech.string_for_config())
438 else:
439 d['str_source'] = '0'
441 template = '''# autogenerated QSEISR input by qseis2d.py
442# This is the input file of FORTRAN77 program "QseisR" for calculation of
443# synthetic seismograms using the given f-k spectra of incident seismic
444# waves at the receiver-site basement, where the f-k spectra should be
445# prepared by the "QseisS" code.
446#
447# by
448# Rongjiang Wang <wang@gfz-potsdam.de>
449# GeoForschungsZentrum Potsdam
450# Telegrafenberg, D-14473 Potsdam, Germany
451#
452# Modified from qseis2006, Potsdam, Dec. 2014
453#
454#------------------------------------------------------------------------------
455# SOURCE PARAMETERS
456# =================
457# 1. epicenter (lat[deg], lon[deg])
458# 2. moment tensor in N*m: Mxx, Myy, Mzz, Mxy, Myz, Mzx
459# Note: x is northward, y is eastward and z is downard
460# conversion from CMT system:
461# Mxx = Mtt, Myy= Mpp, Mzz = Mrr, Mxy = -Mtp, Myz = -Mrp, Mzx = Mrt
462# 3. file of f-k spectra of incident waves
463# 4. source duration [s], and selection of source time functions, i.e.,
464# source wavelet (0 = user's own wavelet; 1 = default wavelet: normalized
465# square half-sinusoid)
466# 5. if user's own wavelet is selected, then number of the wavelet time samples
467# (<= 1024), followed by
468# 6. the equidistant wavelet time series (no comment lines between the time
469# series)
470#------------------------------------------------------------------------------
471 %(str_source_location)s |dble(2);
472 %(str_source)s |dble(6);
473 '%(fk_path)s' |char;
474 %(wavelet_duration)e %(wavelet_type)i %(str_w_samples)s |dble, int, dbls;
475#------------------------------------------------------------------------------
476# RECEIVER PARAMETERS
477# ===================
478# 1. station location (lat[deg], lon[deg], depth[km])
479# Note: the epicentral distance should not exceed the max. distance used
480# for generating the f-k spectra
481# 2. time reduction[s]
482# 3. order of bandpass filter(if <= 0, then no filter will be used), lower
483# and upper cutoff frequences[Hz]
484# 4. selection of output format (1 = Down/Radial/Azimuthal, 2 = East/North/Up)
485# 5. output file of velocity seismograms
486# 6. number of datalines representing the layered receiver-site structure, and
487# selection of surface condition (0 = free surface, 1 = without free
488# surface, i.e., upper halfspace is replaced by medium of the 1. layer)
489# 7. ... layered structure model
490#------------------------------------------------------------------------------
491 %(str_receiver)s |dble(3);
492 %(time_reduction)e |dble;
493 %(str_receiver_filter)s |int, dble(2);
494 %(output_format)i |int;
495 %(str_output_filename)s |char;
496#------------------------------------------------------------------------------
497 %(n_model_receiver_lines)i |int: no_model_lines;
498#------------------------------------------------------------------------------
499# MULTILAYERED MODEL PARAMETERS (shallow receiver-site structure)
500# ===============================================================
501# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs
502#------------------------------------------------------------------------------
503%(model_receiver_lines)s
504#-----------------------END OF INPUT PARAMETERS--------------------------------
505#
506#-Requirements to use QSEIS2d: ------------------------------------------------
507# (1) Teleseismic body waves with penetration depth much larger than the
508# receiver-site basement depth
509# (2) The last layer parameters of the receiver-site structure should be
510# identical with that of the source-site model at the depth which is
511# defined as the common basement depth
512# (3) The cutoff frequency should be high enough for separating different
513# wave types.
514''' # noqa
515 return (template % d).encode('ascii')
518class QSeis2dConfig(Object):
519 '''
520 Combined config object for QSeisS and QSeisR.
522 This config object should contain all settings which cannot be derived from
523 the backend-independant Pyrocko GF Store config.
524 '''
526 qseis_s_config = QSeisSConfig.T(default=QSeisSConfig.D())
527 qseis_r_config = QSeisRConfig.T(default=QSeisRConfig.D())
528 qseis2d_version = String.T(default='2014')
530 time_region = Tuple.T(2, Timing.T(), default=default_time_region)
531 cut = Tuple.T(2, Timing.T(), optional=True)
532 fade = Tuple.T(4, Timing.T(), optional=True)
533 relevel_with_fade_in = Bool.T(default=False)
535 gf_directory = String.T(default='qseis2d_green')
538class QSeis2dError(gf.store.StoreError):
539 pass
542class Interrupted(gf.store.StoreError):
543 def __str__(self):
544 return 'Interrupted.'
547class QSeisSRunner(object):
548 '''
549 Takes QSeis2dConfigFull or QSeisSConfigFull objects, runs the program.
550 '''
551 def __init__(self, tmp, keep_tmp=False):
552 self.tempdir = mkdtemp(prefix='qseisSrun-', dir=tmp)
553 self.keep_tmp = keep_tmp
554 self.config = None
556 def run(self, config):
558 if isinstance(config, QSeis2dConfig):
559 config = QSeisSConfig(**config.qseis_s_config)
561 self.config = config
563 input_fn = op.join(self.tempdir, 'input')
565 with open(input_fn, 'wb') as f:
566 input_str = config.string_for_config()
567 logger.debug('===== begin qseisS input =====\n'
568 '%s===== end qseisS input =====' % input_str.decode())
569 f.write(input_str)
571 program = program_bins['qseis2d.qseisS%s' % config.qseiss_version]
573 old_wd = os.getcwd()
574 os.chdir(self.tempdir)
576 interrupted = []
578 def signal_handler(signum, frame):
579 os.kill(proc.pid, signal.SIGTERM)
580 interrupted.append(True)
582 original = signal.signal(signal.SIGINT, signal_handler)
583 try:
584 try:
585 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
587 except OSError:
588 os.chdir(old_wd)
589 raise QSeis2dError('could not start qseisS: "%s"' % program)
591 (output_str, error_str) = proc.communicate(b'input\n')
593 finally:
594 signal.signal(signal.SIGINT, original)
596 if interrupted:
597 raise KeyboardInterrupt()
599 logger.debug('===== begin qseisS output =====\n'
600 '%s===== end qseisS output =====' % output_str.decode())
602 errmess = []
603 if proc.returncode != 0:
604 errmess.append(
605 'qseisS had a non-zero exit state: %i' % proc.returncode)
607 if error_str:
608 errmess.append('qseisS emitted something via stderr')
610 if output_str.lower().find(b'error') != -1:
611 errmess.append("the string 'error' appeared in qseisS output")
613 if errmess:
614 self.keep_tmp = True
616 os.chdir(old_wd)
617 raise QSeis2dError('''
618===== begin qseisS input =====
619%s===== end qseisS input =====
620===== begin qseisS output =====
621%s===== end qseisS output =====
622===== begin qseisS error =====
623%s===== end qseisS error =====
624%s
625qseisS has been invoked as "%s"
626in the directory %s'''.lstrip() % (
627 input_str.decode(),
628 output_str.decode(),
629 error_str.decode(),
630 '\n'.join(errmess), program,
631 self.tempdir))
633 self.qseiss_output = output_str
634 self.qseiss_error = error_str
636 os.chdir(old_wd)
638 def __del__(self):
639 if self.tempdir:
640 if not self.keep_tmp:
641 shutil.rmtree(self.tempdir)
642 self.tempdir = None
643 else:
644 logger.warn(
645 'not removing temporary directory: %s' % self.tempdir)
648class QSeisRRunner(object):
649 '''
650 Takes QSeis2dConfig or QSeisRConfigFull objects, runs the program and
651 reads the output.
652 '''
653 def __init__(self, tmp=None, keep_tmp=False):
654 self.tempdir = mkdtemp(prefix='qseisRrun-', dir=tmp)
655 self.keep_tmp = keep_tmp
656 self.config = None
658 def run(self, config):
659 if isinstance(config, QSeis2dConfig):
660 config = QSeisRConfigFull(**config.qseis_r_config)
662 self.config = config
664 input_fn = op.join(self.tempdir, 'input')
666 with open(input_fn, 'wb') as f:
667 input_str = config.string_for_config()
668 old_wd = os.getcwd()
669 os.chdir(self.tempdir)
670 logger.debug('===== begin qseisR input =====\n'
671 '%s===== end qseisR input =====' % input_str.decode())
672 f.write(input_str)
674 program = program_bins['qseis2d.qseisR%s' % config.qseisr_version]
676 interrupted = []
678 def signal_handler(signum, frame):
679 os.kill(proc.pid, signal.SIGTERM)
680 interrupted.append(True)
682 original = signal.signal(signal.SIGINT, signal_handler)
683 try:
684 try:
685 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
687 except OSError:
688 os.chdir(old_wd)
689 raise QSeis2dError('could not start qseisR: "%s"' % program)
691 (output_str, error_str) = proc.communicate(b'input\n')
693 finally:
694 signal.signal(signal.SIGINT, original)
696 if interrupted:
697 raise KeyboardInterrupt()
699 logger.debug('===== begin qseisR output =====\n'
700 '%s===== end qseisR output =====' % output_str.decode())
702 errmess = []
703 if proc.returncode != 0:
704 errmess.append(
705 'qseisR had a non-zero exit state: %i' % proc.returncode)
707 if error_str:
708 errmess.append('qseisR emitted something via stderr')
710 if output_str.lower().find(b'error') != -1:
711 errmess.append("the string 'error' appeared in qseisR output")
713 if errmess:
714 self.keep_tmp = True
716 os.chdir(old_wd)
717 raise QSeis2dError('''
718===== begin qseisR input =====
719%s===== end qseisR input =====
720===== begin qseisR output =====
721%s===== end qseisR output =====
722===== begin qseisR error =====
723%s===== end qseisR error =====
724%s
725qseisR has been invoked as "%s"
726in the directory %s'''.lstrip() % (
727 input_str.decode(),
728 output_str.decode(),
729 error_str.decode(),
730 '\n'.join(errmess), program,
731 self.tempdir))
733 self.qseisr_output = output_str
734 self.qseisr_error = error_str
736 os.chdir(old_wd)
738 def get_traces(self):
740 fn = self.config.get_output_filename(self.tempdir)
741 data = num.loadtxt(fn, skiprows=1, dtype=float)
742 nsamples, ntraces = data.shape
743 deltat = (data[-1, 0] - data[0, 0]) / (nsamples - 1)
744 toffset = data[0, 0]
746 tred = self.config.time_reduction
747 rec = self.config.receiver
748 tmin = rec.tstart + toffset + deltat + tred
750 traces = []
751 for itrace, comp in enumerate(self.config.components):
752 # qseis2d gives velocity-integrate to displacement
753 # integration removes one sample, add it again in front
754 displ = cumtrapz(num.concatenate(
755 (num.zeros(1), data[:, itrace + 1])), dx=deltat)
757 tr = trace.Trace(
758 '', '%04i' % itrace, '', comp,
759 tmin=tmin, deltat=deltat, ydata=displ,
760 meta=dict(distance=rec.distance,
761 azimuth=0.0))
763 traces.append(tr)
765 return traces
767 def __del__(self):
768 if self.tempdir:
769 if not self.keep_tmp:
770 shutil.rmtree(self.tempdir)
771 self.tempdir = None
772 else:
773 logger.warn(
774 'not removing temporary directory: %s' % self.tempdir)
777class QSeis2dGFBuilder(gf.builder.Builder):
778 nsteps = 2
780 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
781 force=False):
783 self.store = gf.store.Store(store_dir, 'w')
785 storeconf = self.store.config
787 if step == 0:
788 block_size = (1, 1, storeconf.ndistances)
789 else:
790 if block_size is None:
791 block_size = (1, 1, 1) # QSeisR does only allow one receiver
793 if len(storeconf.ns) == 2:
794 block_size = block_size[1:]
796 gf.builder.Builder.__init__(
797 self, storeconf, step, block_size=block_size)
799 baseconf = self.store.get_extra('qseis2d')
801 conf_s = QSeisSConfigFull(**baseconf.qseis_s_config.items())
802 conf_r = QSeisRConfigFull(**baseconf.qseis_r_config.items())
804 conf_s.earthmodel_1d = storeconf.earthmodel_1d
805 if storeconf.earthmodel_receiver_1d is not None:
806 conf_r.earthmodel_receiver_1d = \
807 storeconf.earthmodel_receiver_1d
809 else:
810 conf_r.earthmodel_receiver_1d = \
811 storeconf.earthmodel_1d.extract(
812 depth_max='moho')
813 # depth_max=conf_s.receiver_basement_depth*km)
815 deltat = 1.0 / self.gf_config.sample_rate
817 if 'time_window_min' not in shared:
818 d = self.store.make_timing_params(
819 baseconf.time_region[0], baseconf.time_region[1],
820 force=force)
822 shared['time_window_min'] = float(
823 num.ceil(d['tlenmax'] / self.gf_config.sample_rate) *
824 self.gf_config.sample_rate)
825 shared['time_reduction'] = d['tmin_vred']
827 time_window_min = shared['time_window_min']
829 conf_s.nsamples = nextpow2(int(round(time_window_min / deltat)) + 1)
830 conf_s.time_window = (conf_s.nsamples - 1) * deltat
831 conf_r.time_reduction = shared['time_reduction']
833 if step == 0:
834 if 'slowness_window' not in shared:
835 if conf_s.calc_slowness_window:
836 phases = [
837 storeconf.tabulated_phases[i].phases
838 for i in range(len(
839 storeconf.tabulated_phases))]
841 all_phases = []
842 map(all_phases.extend, phases)
844 mean_source_depth = num.mean((
845 storeconf.source_depth_min,
846 storeconf.source_depth_max))
848 arrivals = conf_s.earthmodel_1d.arrivals(
849 phases=all_phases,
850 distances=num.linspace(
851 conf_s.receiver_min_distance,
852 conf_s.receiver_max_distance,
853 100) * cake.m2d,
854 zstart=mean_source_depth)
856 ps = num.array(
857 [arrivals[i].p for i in range(len(arrivals))])
859 slownesses = ps / (cake.r2d * cake.d2m / km)
861 shared['slowness_window'] = (0.,
862 0.,
863 1.1 * float(slownesses.max()),
864 1.3 * float(slownesses.max()))
866 else:
867 shared['slowness_window'] = conf_s.slowness_window
869 conf_s.slowness_window = shared['slowness_window']
871 self.qseis_s_config = conf_s
872 self.qseis_r_config = conf_r
873 self.qseis_baseconf = baseconf
875 self.tmp = tmp
876 if self.tmp is not None:
877 util.ensuredir(self.tmp)
879 util.ensuredir(baseconf.gf_directory)
881 def cleanup(self):
882 self.store.close()
884 def work_block(self, iblock):
885 if len(self.store.config.ns) == 2:
886 (sz, firstx), (sz, lastx), (ns, nx) = \
887 self.get_block_extents(iblock)
889 rz = self.store.config.receiver_depth
890 else:
891 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
892 self.get_block_extents(iblock)
894 source_depth = float(sz / km)
895 conf_s = copy.deepcopy(self.qseis_s_config)
896 conf_r = copy.deepcopy(self.qseis_r_config)
898 gf_directory = op.abspath(self.qseis_baseconf.gf_directory)
900 fk_path = op.join(gf_directory, 'green_%.3fkm.fk' % source_depth)
901 info_path = op.join(gf_directory, 'green_%.3fkm.info' % source_depth)
903 conf_s.fk_path = fk_path
904 conf_s.info_path = info_path
906 conf_r.fk_path = fk_path
907 conf_r.info_path = info_path
909 if self.step == 0 and os.path.isfile(fk_path):
910 logger.info('Skipping step %i / %i, block %i / %i'
911 '(GF already exists)' %
912 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
913 return
915 logger.info(
916 'Starting step %i / %i, block %i / %i' %
917 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
919 dx = self.gf_config.distance_delta
920 conf_r.wavelet_duration = 0.001 * self.gf_config.sample_rate
922 if self.step == 0:
923 conf_s.source_depth = source_depth
924 runner = QSeisSRunner(tmp=self.tmp)
925 runner.run(conf_s)
927 else:
928 conf_r.receiver = QSeisRReceiver(lat=90 - firstx * cake.m2d,
929 lon=180.,
930 tstart=0.0,
931 distance=firstx)
932 conf_r.source = QSeis2dSource(lat=90 - 0.001 * dx * cake.m2d,
933 lon=0.0,
934 depth=source_depth)
936 runner = QSeisRRunner(tmp=self.tmp)
938 mmt1 = (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
939 {'r': (0, 1), 't': (3, 1), 'z': (5, 1)})
940 mmt2 = (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
941 {'r': (1, 1), 't': (4, 1), 'z': (6, 1)})
942 mmt3 = (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
943 {'r': (2, 1), 'z': (7, 1)})
944 mmt4 = (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
945 {'r': (8, 1), 'z': (9, 1)})
947 gfmapping = [mmt1, mmt2, mmt3, mmt4]
949 for mt, gfmap in gfmapping:
950 if mt:
951 m = mt.m()
952 f = float
953 conf_r.source_mech = QSeisRSourceMechMT(
954 mnn=f(m[0, 0]), mee=f(m[1, 1]), mdd=f(m[2, 2]),
955 mne=f(m[0, 1]), mnd=f(m[0, 2]), med=f(m[1, 2]))
956 else:
957 conf_r.source_mech = None
959 if conf_r.source_mech is not None:
960 runner.run(conf_r)
962 rawtraces = runner.get_traces()
964 interrupted = []
966 def signal_handler(signum, frame):
967 interrupted.append(True)
969 original = signal.signal(signal.SIGINT, signal_handler)
970 self.store.lock()
971 duplicate_inserts = 0
972 try:
973 for itr, tr in enumerate(rawtraces):
974 if tr.channel not in gfmap:
975 continue
977 x = tr.meta['distance']
978 if x > firstx + (nx - 1) * dx:
979 continue
981 ig, factor = gfmap[tr.channel]
983 if len(self.store.config.ns) == 2:
984 args = (sz, x, ig)
985 else:
986 args = (rz, sz, x, ig)
988 if self.qseis_baseconf.cut:
989 tmin = self.store.t(
990 self.qseis_baseconf.cut[0], args[:-1])
991 tmax = self.store.t(
992 self.qseis_baseconf.cut[1], args[:-1])
994 if None in (tmin, tmax):
995 continue
997 tr.chop(tmin, tmax)
999 tmin = tr.tmin
1000 tmax = tr.tmax
1002 if self.qseis_baseconf.fade:
1003 ta, tb, tc, td = [
1004 self.store.t(v, args[:-1])
1005 for v in self.qseis_baseconf.fade]
1007 if None in (ta, tb, tc, td):
1008 continue
1010 if not (ta <= tb and tb <= tc and tc <= td):
1011 raise QSeis2dError(
1012 'invalid fade configuration')
1014 t = tr.get_xdata()
1015 fin = num.interp(t, [ta, tb], [0., 1.])
1016 fout = num.interp(t, [tc, td], [1., 0.])
1017 anti_fin = 1. - fin
1018 anti_fout = 1. - fout
1020 y = tr.ydata
1022 sum_anti_fin = num.sum(anti_fin)
1023 sum_anti_fout = num.sum(anti_fout)
1025 if sum_anti_fin != 0.0:
1026 yin = num.sum(anti_fin * y) / sum_anti_fin
1027 else:
1028 yin = 0.0
1030 if sum_anti_fout != 0.0:
1031 yout = num.sum(anti_fout * y) / sum_anti_fout
1032 else:
1033 yout = 0.0
1035 y2 = anti_fin * yin + \
1036 fin * fout * y + \
1037 anti_fout * yout
1039 if self.qseis_baseconf.relevel_with_fade_in:
1040 y2 -= yin
1042 tr.set_ydata(y2)
1044 gf_tr = gf.store.GFTrace.from_trace(tr)
1045 gf_tr.data *= factor
1047 try:
1048 self.store.put(args, gf_tr)
1049 except gf.store.DuplicateInsert:
1050 duplicate_inserts += 1
1052 finally:
1053 if duplicate_inserts:
1054 logger.warn('%i insertions skipped (duplicates)' %
1055 duplicate_inserts)
1057 self.store.unlock()
1058 signal.signal(signal.SIGINT, original)
1060 if interrupted:
1061 raise KeyboardInterrupt()
1063 logger.info(
1064 'Done with step %i / %i, block %i / %i' %
1065 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
1068def init(store_dir, variant):
1069 if variant is None:
1070 variant = '2014'
1072 if variant not in ('2014'):
1073 raise gf.store.StoreError('unsupported variant: %s' % variant)
1075 modelling_code_id = 'qseis2d.%s' % variant
1077 qseis2d = QSeis2dConfig()
1079 qseis2d.time_region = (
1080 gf.meta.Timing('begin-50'),
1081 gf.meta.Timing('end+100'))
1083 qseis2d.cut = (
1084 gf.meta.Timing('begin-50'),
1085 gf.meta.Timing('end+100'))
1087 qseis2d.qseis_s_config.sw_flat_earth_transform = 1
1089 store_id = os.path.basename(os.path.realpath(store_dir))
1091 config = gf.meta.ConfigTypeA(
1093 id=store_id,
1094 ncomponents=10,
1095 sample_rate=0.2,
1096 receiver_depth=0 * km,
1097 source_depth_min=10 * km,
1098 source_depth_max=20 * km,
1099 source_depth_delta=10 * km,
1100 distance_min=100 * km,
1101 distance_max=1000 * km,
1102 distance_delta=10 * km,
1103 earthmodel_1d=cake.load_model().extract(depth_max='cmb'),
1104 modelling_code_id=modelling_code_id,
1105 tabulated_phases=[
1106 gf.meta.TPDef(
1107 id='begin',
1108 definition='p,P,p\\,P\\,Pv_(cmb)p'),
1109 gf.meta.TPDef(
1110 id='end',
1111 definition='2.5'),
1112 gf.meta.TPDef(
1113 id='P',
1114 definition='!P'),
1115 gf.meta.TPDef(
1116 id='S',
1117 definition='!S'),
1118 gf.meta.TPDef(
1119 id='p',
1120 definition='!p'),
1121 gf.meta.TPDef(
1122 id='s',
1123 definition='!s')])
1125 config.validate()
1126 return gf.store.Store.create_editables(
1127 store_dir, config=config, extra={'qseis2d': qseis2d})
1130def build(store_dir, force=False, nworkers=None, continue_=False, step=None,
1131 iblock=None):
1133 return QSeis2dGFBuilder.build(
1134 store_dir, force=force, nworkers=nworkers, continue_=continue_,
1135 step=step, iblock=iblock)