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
16import os.path as op
17from scipy.integrate import cumtrapz
19from pyrocko.moment_tensor import MomentTensor, symmat6
20from pyrocko.guts import Float, Int, Tuple, List, Bool, Object, String
21from pyrocko import trace, util, cake, gf
23km = 1e3
25guts_prefix = 'pf'
27Timing = gf.meta.Timing
30logger = logging.getLogger('pyrocko.fomosto.qseis2d')
32# how to call the programs
33program_bins = {
34 'qseis2d.qseisS2014': 'fomosto_qseisS2014',
35 'qseis2d.qseisR2014': 'fomosto_qseisR2014',
36}
38qseis2d_components = {
39 1: 'z r t'.split(),
40 2: 'e n u'.split(),
41}
43# defaults
44default_gf_directory = 'qseis2d_green'
45default_fk_basefilename = 'green'
46default_source_depth = 10.0
47default_time_region = (Timing('-10'), Timing('+890'))
48default_slowness_window = (0.0, 0.0, 0.2, 0.25)
51def have_backend():
52 for cmd in [[exe] for exe in program_bins.values()]:
53 try:
54 p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE)
55 (stdout, stderr) = p.communicate()
57 except OSError:
58 return False
60 return True
63def nextpow2(i):
64 return 2 ** int(math.ceil(math.log(i) / math.log(2.)))
67def str_float_vals(vals):
68 return ' '.join('%e' % val for val in vals)
71def str_int_vals(vals):
72 return ' '.join('%i' % val for val in vals)
75def str_str_vals(vals):
76 return ' '.join("'%s'" % val for val in vals)
79def scl(cs):
80 if not cs:
81 return '\n#'
83 return '\n' + ' '.join('(%e,%e)' % (c.real, c.imag) for c in cs)
86def cake_model_to_config(mod):
87 k = 1000.
88 srows = []
89 for i, row in enumerate(mod.to_scanlines()):
90 depth, vp, vs, rho, qp, qs = row
91 row = [depth / k, vp / k, vs / k, rho / k, qp, qs]
92 srows.append('%i %s' % (i + 1, str_float_vals(row)))
94 return '\n'.join(srows), len(srows)
97class QSeis2dSource(Object):
98 lat = Float.T(default=10.)
99 lon = Float.T(default=0.)
100 depth = Float.T(default=default_source_depth)
102 def string_for_config(self):
103 return '%(lat)e %(lon)15e ' % self.__dict__
106class QSeisRSourceMech(Object):
107 pass
110class QSeisRSourceMechMT(QSeisRSourceMech):
111 mnn = Float.T(default=1.0)
112 mee = Float.T(default=1.0)
113 mdd = Float.T(default=1.0)
114 mne = Float.T(default=0.0)
115 mnd = Float.T(default=0.0)
116 med = Float.T(default=0.0)
118 def string_for_config(self):
119 return '%(mnn)e %(mee)15e %(mdd)15e ' \
120 '%(mne)15e %(med)15e %(mnd)15e ' % self.__dict__
123class QSeisSPropagationFilter(Object):
124 min_depth = Float.T(default=0.0)
125 max_depth = Float.T(default=0.0)
126 filtered_phase = Int.T(default=0)
128 def string_for_config(self):
129 return '%(min_depth)15e %(max_depth)15e ' \
130 '%(filtered_phase)i' % self.__dict__
133class QSeisRBandpassFilter(Object):
134 order = Int.T(default=-1)
135 lower_cutoff = Float.T(default=0.1) # [Hz]
136 upper_cutoff = Float.T(default=10.0)
138 def string_for_config(self):
139 return '%(order)i %(lower_cutoff)5e %(upper_cutoff)5e' % self.__dict__
142class QSeisSConfig(Object):
144 qseiss_version = String.T(default='2014')
146 calc_slowness_window = Int.T(default=1)
147 slowness_window = Tuple.T(4, optional=True)
148 wavenumber_sampling = Float.T(default=2.5)
149 aliasing_suppression_factor = Float.T(default=0.01)
151 filter_shallow_paths = Int.T(default=0)
152 filter_shallow_paths_depth = Float.T(default=0.0)
153 propagation_filters = List.T(QSeisSPropagationFilter.T())
155 sw_flat_earth_transform = Int.T(default=0)
157 gradient_resolution_vp = Float.T(default=0.0)
158 gradient_resolution_vs = Float.T(default=0.0)
159 gradient_resolution_density = Float.T(default=0.0)
161 def items(self):
162 return dict(self.T.inamevals(self))
165class QSeisSConfigFull(QSeisSConfig):
167 time_window = Float.T(default=900.0)
169 source_depth = Float.T(default=10.0)
171 receiver_basement_depth = Float.T(default=35.0) # [km]
172 receiver_min_distance = Float.T(default=1000.0) # [km]
173 receiver_max_distance = Float.T(default=10000.0) # [km]
174 nsamples = Int.T(default=256)
176 info_path = String.T(default='green.info')
177 fk_path = String.T(default='green.fk')
179 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True)
181 @staticmethod
182 def example():
183 conf = QSeisSConfigFull()
184 conf.source_depth = 15.
185 conf.receiver_basement_depth = 35.
186 conf.receiver_max_distance = 2000.
187 conf.earthmodel_1d = cake.load_model().extract(depth_max='cmb')
188 conf.sw_flat_earth_transform = 1
189 return conf
191 def string_for_config(self):
193 def aggregate(xx):
194 return len(xx), '\n'.join(
195 [''] + [x.string_for_config() for x in xx])
197 assert self.earthmodel_1d is not None
198 assert self.slowness_window is not None or self.calc_slowness_window,\
199 "'slowness window' undefined and 'calc_slowness_window'=0"
201 d = self.__dict__.copy()
203 model_str, nlines = cake_model_to_config(self.earthmodel_1d)
204 d['n_model_lines'] = nlines
205 d['model_lines'] = model_str
207 if not self.slowness_window:
208 d['str_slowness_window'] = str_float_vals(default_slowness_window)
209 else:
210 d['str_slowness_window'] = str_float_vals(self.slowness_window)
212 d['n_depth_ranges'], d['str_depth_ranges'] = \
213 aggregate(self.propagation_filters)
215 template = '''# autogenerated QSEISS input by qseis2d.py
216# This is the input file of FORTRAN77 program "QseisS" for calculation of
217# f-k spectra of upgoing seismic waves at the reveiver-site basement.
218#
219# by
220# Rongjiang Wang <wang@gfz-potsdam.de>
221# GeoForschungsZentrum Potsdam
222# Telegrafenberg, D-14473 Potsdam, Germany
223#
224# Modified from qseis2006, Potsdam, Dec. 2014
225#
226# SOURCE DEPTH
227# ============
228# 1. source depth [km]
229#------------------------------------------------------------------------------
230 %(source_depth)e |dble;
231#------------------------------------------------------------------------------
232#
233# RECEIVER-SITE PARAMETERS
234# ========================
235# 1. receiver-site basement depth [km]
236# 2. max. epicental distance [km]
237#------------------------------------------------------------------------------
238 %(receiver_basement_depth)e |dble;
239 %(receiver_max_distance)e |dble;
240#------------------------------------------------------------------------------
241# TIME SAMPLING PARAMETERS
242# ========================
243# 1. length of time window [sec]
244# 2. number of time samples (<= 2*nfmax in qsglobal.h)
245#------------------------------------------------------------------------------
246 %(time_window)e |dble: t_window;
247 %(nsamples)i |int: no_t_samples;
248#------------------------------------------------------------------------------
249# SLOWNESS WINDOW PARAMETERS
250# ==========================
251# 1. the low and high slowness cut-offs [s/km] with tapering: 0 < slw1 < slw2
252# defining the bounds of cosine taper at the lower end, and 0 < slw3 < slw4
253# defining the bounds of cosine taper at the higher end.
254# 2. slowness sampling rate (1.0 = sampling with the Nyquist slowness, 2.0 =
255# sampling with twice higher than the Nyquist slowness, and so on: the
256# larger this parameter, the smaller the space-domain aliasing effect, but
257# also the more computation effort);
258# 3. the factor for suppressing time domain aliasing (> 0 and <= 1).
259#------------------------------------------------------------------------------
260 %(str_slowness_window)s |dble: slw(1-4);
261 %(wavenumber_sampling)e |dble: sample_rate;
262 %(aliasing_suppression_factor)e |dble: supp_factor;
263#------------------------------------------------------------------------------
264# OPTIONS FOR PARTIAL SOLUTIONS
265# =============================
266# 1. switch for filtering waves with a shallow penetration depth (concerning
267# their whole trace from source to receiver), penetration depth limit [km]
268# (Note: if this option is selected, waves whose travel path never exceeds
269# the given depth limit will be filtered ("seismic nuting"). the condition
270# for selecting this filter is that the given shallow path depth limit
271# should be larger than both source and receiver depth.)
272# 2. number of depth ranges where the following selected up/down-sp2oing P or
273# SV waves should be filtered
274# 3. the 1. depth range: upper and lower depth [km], switch for filtering P
275# or SV wave in this depth range:
276# switch no: 1 2 3 4 other
277# filtered phase: P(up) P(down) SV(up) SV(down) Error
278# 4. the 2. ...
279# (Note: the partial solution options are useful tools to increase the
280# numerical resolution of desired wave phases. especially when the desired
281# phases are much smaller than the undesired phases, these options should
282# be selected and carefully combined.)
283#------------------------------------------------------------------------------
284 %(filter_shallow_paths)i %(filter_shallow_paths_depth)e
285 %(n_depth_ranges)i %(str_depth_ranges)s |int, dble, dble, int;
286#------------------------------------------------------------------------------
287# OUTPUT FILE FOR F-K SPECTRA of GREEN'S FUNCTIONS
288# ================================================
289# 1. info-file of Green's functions (ascii including f-k sampling parameters)
290# 2. file name of Green's functions (binary files including explosion, strike
291# -slip, dip-slip and clvd sources)
292#------------------------------------------------------------------------------
293 '%(info_path)s'
294 '%(fk_path)s'
295#------------------------------------------------------------------------------
296# GLOBAL MODEL PARAMETERS
297# =======================
298# 1. switch for flat-earth-transform
299# 2. gradient resolution [percent] of vp, vs, and rho (density), if <= 0, then
300# default values (depending on wave length at cut-off frequency) will be
301# used
302#------------------------------------------------------------------------------
303 %(sw_flat_earth_transform)i |int: sw_flat_earth_transform;
304 %(gradient_resolution_vp)e %(gradient_resolution_vs)e %(gradient_resolution_density)e
305#------------------------------------------------------------------------------
306# SOURCE-SITE LAYERED EARTH MODEL
307# ===============================
308# 1. number of data lines of the layered model
309#------------------------------------------------------------------------------
310 %(n_model_lines)i |int: no_model_lines;
311#------------------------------------------------------------------------------
312# no depth[km] vp[km/s] vs[km/s] rho[g/cm^3] qp qs
313#------------------------------------------------------------------------------
314%(model_lines)s ;
315#-----------------------END OF INPUT PARAMETERS--------------------------------
317Glossary
319SLOWNESS: The slowness is the inverse of apparent wave velocity = sin(i)/v with
320i = incident angle and v = true wave velocity.
322SUPPRESSION OF TIME-DOMAIN ALIASING: The suppression of the time domain aliasing
323is achieved by using the complex frequency technique. The suppression factor
324should be a value between 0 and 1. If this factor is set to 0.1, for example, the
325aliasing phase at the reduced time begin is suppressed to 10 percent.
327MODEL PARAMETER GRADIENT RESOLUTION: Layers with a constant gradient will be
328discretized with a number of homogeneous sublayers. The gradient resolutions are
329then used to determine the maximum allowed thickness of the sublayers. If the
330resolutions of Vp, Vs and Rho (density) require different thicknesses, the
331smallest is first chosen. If this is even smaller than 1 percent of the
332characteristic wavelength, then the latter is taken finally for the sublayer
333thickness.
334''' # noqa
335 return (template % d).encode('ascii')
338class QSeisRReceiver(Object):
339 lat = Float.T(default=10.0)
340 lon = Float.T(default=0.0)
341 depth = Float.T(default=0.0)
342 tstart = Float.T(default=0.0)
343 distance = Float.T(default=0.0)
345 def string_for_config(self):
346 return '%(lat)e %(lon)15e %(depth)15e' % self.__dict__
349class QSeisRConfig(Object):
351 qseisr_version = String.T(default='2014')
353 receiver_filter = QSeisRBandpassFilter.T(optional=True)
355 wavelet_duration = Float.T(default=0.001) # [s]
356 wavelet_type = Int.T(default=1)
357 user_wavelet_samples = List.T(Float.T())
359 def items(self):
360 return dict(self.T.inamevals(self))
363class QSeisRConfigFull(QSeisRConfig):
365 source = QSeis2dSource(depth=default_source_depth) # [lat, lon(deg)]
366 receiver = QSeisRReceiver() # [lat, lon(deg)]
368 source_mech = QSeisRSourceMech.T(
369 optional=True,
370 default=QSeisRSourceMechMT.D())
372 time_reduction = Float.T(default=0.0)
374 info_path = String.T(default='green.info')
375 fk_path = String.T(default='green.fk')
377 output_format = Int.T(default=1) # 1/2 components in [Z, R, T]/[E, N, U]
378 output_filename = String.T(default='seis.dat')
380 earthmodel_receiver_1d = gf.meta.Earthmodel1D.T(optional=True)
382 @staticmethod
383 def example():
384 conf = QSeisRConfigFull()
385 conf.source = QSeis2dSource(lat=-80.5, lon=90.1)
386 conf.receiver_location = QSeisRReceiver(lat=13.4, lon=240.5, depth=0.0)
387 conf.time_reduction = 10.0
388 conf.earthmodel_receiver_1d = cake.load_model().extract(
389 depth_max='moho')
390 return conf
392 @property
393 def components(self):
394 return qseis2d_components[self.output_format]
396 def get_output_filename(self, rundir):
397 return op.join(rundir, self.output_filename)
399 def string_for_config(self):
401 def aggregate(xx):
402 return len(xx), '\n'.join(
403 [''] + [x.string_for_config() for x in xx])
405 assert self.earthmodel_receiver_1d is not None
407 d = self.__dict__.copy()
409# Actually not existing anymore in code
410# d['sw_surface'] = 0 # 0-free-surface, 1-no fre surface
412 d['str_source_location'] = self.source.string_for_config()
414 d['str_receiver'] = self.receiver.string_for_config()
416 d['str_output_filename'] = "'%s'" % self.output_filename
418 model_str, nlines = cake_model_to_config(self.earthmodel_receiver_1d)
420 d['n_model_receiver_lines'] = nlines
421 d['model_receiver_lines'] = model_str
423 if self.wavelet_type == 0: # user wavelet
424 d['str_w_samples'] = '\n' \
425 + '%i\n' % len(self.user_wavelet_samples) \
426 + str_float_vals(self.user_wavelet_samples)
427 else:
428 d['str_w_samples'] = ''
430 if self.receiver_filter:
431 d['str_receiver_filter'] = self.receiver_filter.string_for_config()
432 else:
433 d['str_receiver_filter'] = '-1 0.0 200.'
435 if self.source_mech:
436 d['str_source'] = '%s' % (self.source_mech.string_for_config())
437 else:
438 d['str_source'] = '0'
440 template = '''# autogenerated QSEISR input by qseis2d.py
441# This is the input file of FORTRAN77 program "QseisR" for calculation of
442# synthetic seismograms using the given f-k spectra of incident seismic
443# waves at the receiver-site basement, where the f-k spectra should be
444# prepared by the "QseisS" code.
445#
446# by
447# Rongjiang Wang <wang@gfz-potsdam.de>
448# GeoForschungsZentrum Potsdam
449# Telegrafenberg, D-14473 Potsdam, Germany
450#
451# Modified from qseis2006, Potsdam, Dec. 2014
452#
453#------------------------------------------------------------------------------
454# SOURCE PARAMETERS
455# =================
456# 1. epicenter (lat[deg], lon[deg])
457# 2. moment tensor in N*m: Mxx, Myy, Mzz, Mxy, Myz, Mzx
458# Note: x is northward, y is eastward and z is downard
459# conversion from CMT system:
460# Mxx = Mtt, Myy= Mpp, Mzz = Mrr, Mxy = -Mtp, Myz = -Mrp, Mzx = Mrt
461# 3. file of f-k spectra of incident waves
462# 4. source duration [s], and selection of source time functions, i.e.,
463# source wavelet (0 = user's own wavelet; 1 = default wavelet: normalized
464# square half-sinusoid)
465# 5. if user's own wavelet is selected, then number of the wavelet time samples
466# (<= 1024), followed by
467# 6. the equidistant wavelet time series (no comment lines between the time
468# series)
469#------------------------------------------------------------------------------
470 %(str_source_location)s |dble(2);
471 %(str_source)s |dble(6);
472 '%(fk_path)s' |char;
473 %(wavelet_duration)e %(wavelet_type)i %(str_w_samples)s |dble, int, dbls;
474#------------------------------------------------------------------------------
475# RECEIVER PARAMETERS
476# ===================
477# 1. station location (lat[deg], lon[deg], depth[km])
478# Note: the epicentral distance should not exceed the max. distance used
479# for generating the f-k spectra
480# 2. time reduction[s]
481# 3. order of bandpass filter(if <= 0, then no filter will be used), lower
482# and upper cutoff frequences[Hz]
483# 4. selection of output format (1 = Down/Radial/Azimuthal, 2 = East/North/Up)
484# 5. output file of velocity seismograms
485# 6. number of datalines representing the layered receiver-site structure, and
486# selection of surface condition (0 = free surface, 1 = without free
487# surface, i.e., upper halfspace is replaced by medium of the 1. layer)
488# 7. ... layered structure model
489#------------------------------------------------------------------------------
490 %(str_receiver)s |dble(3);
491 %(time_reduction)e |dble;
492 %(str_receiver_filter)s |int, dble(2);
493 %(output_format)i |int;
494 %(str_output_filename)s |char;
495#------------------------------------------------------------------------------
496 %(n_model_receiver_lines)i |int: no_model_lines;
497#------------------------------------------------------------------------------
498# MULTILAYERED MODEL PARAMETERS (shallow receiver-site structure)
499# ===============================================================
500# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs
501#------------------------------------------------------------------------------
502%(model_receiver_lines)s
503#-----------------------END OF INPUT PARAMETERS--------------------------------
504#
505#-Requirements to use QSEIS2d: ------------------------------------------------
506# (1) Teleseismic body waves with penetration depth much larger than the
507# receiver-site basement depth
508# (2) The last layer parameters of the receiver-site structure should be
509# identical with that of the source-site model at the depth which is
510# defined as the common basement depth
511# (3) The cutoff frequency should be high enough for separating different
512# wave types.
513''' # noqa
514 return (template % d).encode('ascii')
517class QSeis2dConfig(Object):
518 '''
519 Combined config object for QSeisS and QSeisR.
521 This config object should contain all settings which cannot be derived from
522 the backend-independant Pyrocko GF Store config.
523 '''
525 qseis_s_config = QSeisSConfig.T(default=QSeisSConfig.D())
526 qseis_r_config = QSeisRConfig.T(default=QSeisRConfig.D())
527 qseis2d_version = String.T(default='2014')
529 time_region = Tuple.T(2, Timing.T(), default=default_time_region)
530 cut = Tuple.T(2, Timing.T(), optional=True)
531 fade = Tuple.T(4, Timing.T(), optional=True)
532 relevel_with_fade_in = Bool.T(default=False)
534 gf_directory = String.T(default='qseis2d_green')
537class QSeis2dError(gf.store.StoreError):
538 pass
541class Interrupted(gf.store.StoreError):
542 def __str__(self):
543 return 'Interrupted.'
546class QSeisSRunner(object):
547 '''
548 Takes QSeis2dConfigFull or QSeisSConfigFull objects, runs the program.
549 '''
550 def __init__(self, tmp, keep_tmp=False):
551 self.tempdir = mkdtemp(prefix='qseisSrun-', dir=tmp)
552 self.keep_tmp = keep_tmp
553 self.config = None
555 def run(self, config):
557 if isinstance(config, QSeis2dConfig):
558 config = QSeisSConfig(**config.qseis_s_config)
560 self.config = config
562 input_fn = op.join(self.tempdir, 'input')
564 with open(input_fn, 'wb') as f:
565 input_str = config.string_for_config()
566 logger.debug('===== begin qseisS input =====\n'
567 '%s===== end qseisS input =====' % input_str.decode())
568 f.write(input_str)
570 program = program_bins['qseis2d.qseisS%s' % config.qseiss_version]
572 old_wd = os.getcwd()
573 os.chdir(self.tempdir)
575 interrupted = []
577 def signal_handler(signum, frame):
578 os.kill(proc.pid, signal.SIGTERM)
579 interrupted.append(True)
581 original = signal.signal(signal.SIGINT, signal_handler)
582 try:
583 try:
584 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
586 except OSError:
587 os.chdir(old_wd)
588 raise QSeis2dError('could not start qseisS: "%s"' % program)
590 (output_str, error_str) = proc.communicate(b'input\n')
592 finally:
593 signal.signal(signal.SIGINT, original)
595 if interrupted:
596 raise KeyboardInterrupt()
598 logger.debug('===== begin qseisS output =====\n'
599 '%s===== end qseisS output =====' % output_str.decode())
601 errmess = []
602 if proc.returncode != 0:
603 errmess.append(
604 'qseisS had a non-zero exit state: %i' % proc.returncode)
606 if error_str:
607 errmess.append('qseisS emitted something via stderr')
609 if output_str.lower().find(b'error') != -1:
610 errmess.append("the string 'error' appeared in qseisS output")
612 if errmess:
613 self.keep_tmp = True
615 os.chdir(old_wd)
616 raise QSeis2dError('''
617===== begin qseisS input =====
618%s===== end qseisS input =====
619===== begin qseisS output =====
620%s===== end qseisS output =====
621===== begin qseisS error =====
622%s===== end qseisS error =====
623%s
624qseisS has been invoked as "%s"
625in the directory %s'''.lstrip() % (
626 input_str.decode(),
627 output_str.decode(),
628 error_str.decode(),
629 '\n'.join(errmess), program,
630 self.tempdir))
632 self.qseiss_output = output_str
633 self.qseiss_error = error_str
635 os.chdir(old_wd)
637 def __del__(self):
638 if self.tempdir:
639 if not self.keep_tmp:
640 shutil.rmtree(self.tempdir)
641 self.tempdir = None
642 else:
643 logger.warning(
644 'not removing temporary directory: %s' % self.tempdir)
647class QSeisRRunner(object):
648 '''
649 Takes QSeis2dConfig or QSeisRConfigFull objects, runs the program and
650 reads the output.
651 '''
652 def __init__(self, tmp=None, keep_tmp=False):
653 self.tempdir = mkdtemp(prefix='qseisRrun-', dir=tmp)
654 self.keep_tmp = keep_tmp
655 self.config = None
657 def run(self, config):
658 if isinstance(config, QSeis2dConfig):
659 config = QSeisRConfigFull(**config.qseis_r_config)
661 self.config = config
663 input_fn = op.join(self.tempdir, 'input')
665 with open(input_fn, 'wb') as f:
666 input_str = config.string_for_config()
667 old_wd = os.getcwd()
668 os.chdir(self.tempdir)
669 logger.debug('===== begin qseisR input =====\n'
670 '%s===== end qseisR input =====' % input_str.decode())
671 f.write(input_str)
673 program = program_bins['qseis2d.qseisR%s' % config.qseisr_version]
675 interrupted = []
677 def signal_handler(signum, frame):
678 os.kill(proc.pid, signal.SIGTERM)
679 interrupted.append(True)
681 original = signal.signal(signal.SIGINT, signal_handler)
682 try:
683 try:
684 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE)
686 except OSError:
687 os.chdir(old_wd)
688 raise QSeis2dError('could not start qseisR: "%s"' % program)
690 (output_str, error_str) = proc.communicate(b'input\n')
692 finally:
693 signal.signal(signal.SIGINT, original)
695 if interrupted:
696 raise KeyboardInterrupt()
698 logger.debug('===== begin qseisR output =====\n'
699 '%s===== end qseisR output =====' % output_str.decode())
701 errmess = []
702 if proc.returncode != 0:
703 errmess.append(
704 'qseisR had a non-zero exit state: %i' % proc.returncode)
706 if error_str:
707 errmess.append('qseisR emitted something via stderr')
709 if output_str.lower().find(b'error') != -1:
710 errmess.append("the string 'error' appeared in qseisR output")
712 if errmess:
713 self.keep_tmp = True
715 os.chdir(old_wd)
716 raise QSeis2dError('''
717===== begin qseisR input =====
718%s===== end qseisR input =====
719===== begin qseisR output =====
720%s===== end qseisR output =====
721===== begin qseisR error =====
722%s===== end qseisR error =====
723%s
724qseisR has been invoked as "%s"
725in the directory %s'''.lstrip() % (
726 input_str.decode(),
727 output_str.decode(),
728 error_str.decode(),
729 '\n'.join(errmess), program,
730 self.tempdir))
732 self.qseisr_output = output_str
733 self.qseisr_error = error_str
735 os.chdir(old_wd)
737 def get_traces(self):
739 fn = self.config.get_output_filename(self.tempdir)
740 data = num.loadtxt(fn, skiprows=1, dtype=float)
741 nsamples, ntraces = data.shape
742 deltat = (data[-1, 0] - data[0, 0]) / (nsamples - 1)
743 toffset = data[0, 0]
745 tred = self.config.time_reduction
746 rec = self.config.receiver
747 tmin = rec.tstart + toffset + deltat + tred
749 traces = []
750 for itrace, comp in enumerate(self.config.components):
751 # qseis2d gives velocity-integrate to displacement
752 # integration removes one sample, add it again in front
753 displ = cumtrapz(num.concatenate(
754 (num.zeros(1), data[:, itrace + 1])), dx=deltat)
756 tr = trace.Trace(
757 '', '%04i' % itrace, '', comp,
758 tmin=tmin, deltat=deltat, ydata=displ,
759 meta=dict(distance=rec.distance,
760 azimuth=0.0))
762 traces.append(tr)
764 return traces
766 def __del__(self):
767 if self.tempdir:
768 if not self.keep_tmp:
769 shutil.rmtree(self.tempdir)
770 self.tempdir = None
771 else:
772 logger.warning(
773 'not removing temporary directory: %s' % self.tempdir)
776class QSeis2dGFBuilder(gf.builder.Builder):
777 nsteps = 2
779 def __init__(self, store_dir, step, shared, block_size=None, tmp=None,
780 force=False):
782 self.store = gf.store.Store(store_dir, 'w')
784 storeconf = self.store.config
786 if step == 0:
787 block_size = (1, 1, storeconf.ndistances)
788 else:
789 if block_size is None:
790 block_size = (1, 1, 1) # QSeisR does only allow one receiver
792 if len(storeconf.ns) == 2:
793 block_size = block_size[1:]
795 gf.builder.Builder.__init__(
796 self, storeconf, step, block_size=block_size)
798 baseconf = self.store.get_extra('qseis2d')
800 conf_s = QSeisSConfigFull(**baseconf.qseis_s_config.items())
801 conf_r = QSeisRConfigFull(**baseconf.qseis_r_config.items())
803 conf_s.earthmodel_1d = storeconf.earthmodel_1d
804 if storeconf.earthmodel_receiver_1d is not None:
805 conf_r.earthmodel_receiver_1d = \
806 storeconf.earthmodel_receiver_1d
808 else:
809 conf_r.earthmodel_receiver_1d = \
810 storeconf.earthmodel_1d.extract(
811 depth_max='moho')
812 # depth_max=conf_s.receiver_basement_depth*km)
814 deltat = 1.0 / self.gf_config.sample_rate
816 if 'time_window_min' not in shared:
817 d = self.store.make_timing_params(
818 baseconf.time_region[0], baseconf.time_region[1],
819 force=force)
821 shared['time_window_min'] = float(
822 num.ceil(d['tlenmax'] / self.gf_config.sample_rate) *
823 self.gf_config.sample_rate)
824 shared['time_reduction'] = d['tmin_vred']
826 time_window_min = shared['time_window_min']
828 conf_s.nsamples = nextpow2(int(round(time_window_min / deltat)) + 1)
829 conf_s.time_window = (conf_s.nsamples - 1) * deltat
830 conf_r.time_reduction = shared['time_reduction']
832 if step == 0:
833 if 'slowness_window' not in shared:
834 if conf_s.calc_slowness_window:
835 phases = [
836 storeconf.tabulated_phases[i].phases
837 for i in range(len(
838 storeconf.tabulated_phases))]
840 all_phases = []
841 map(all_phases.extend, phases)
843 mean_source_depth = num.mean((
844 storeconf.source_depth_min,
845 storeconf.source_depth_max))
847 arrivals = conf_s.earthmodel_1d.arrivals(
848 phases=all_phases,
849 distances=num.linspace(
850 conf_s.receiver_min_distance,
851 conf_s.receiver_max_distance,
852 100) * cake.m2d,
853 zstart=mean_source_depth)
855 ps = num.array(
856 [arrivals[i].p for i in range(len(arrivals))])
858 slownesses = ps / (cake.r2d * cake.d2m / km)
860 shared['slowness_window'] = (0.,
861 0.,
862 1.1 * float(slownesses.max()),
863 1.3 * float(slownesses.max()))
865 else:
866 shared['slowness_window'] = conf_s.slowness_window
868 conf_s.slowness_window = shared['slowness_window']
870 self.qseis_s_config = conf_s
871 self.qseis_r_config = conf_r
872 self.qseis_baseconf = baseconf
874 self.tmp = tmp
875 if self.tmp is not None:
876 util.ensuredir(self.tmp)
878 util.ensuredir(baseconf.gf_directory)
880 def cleanup(self):
881 self.store.close()
883 def work_block(self, iblock):
884 if len(self.store.config.ns) == 2:
885 (sz, firstx), (sz, lastx), (ns, nx) = \
886 self.get_block_extents(iblock)
888 rz = self.store.config.receiver_depth
889 else:
890 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \
891 self.get_block_extents(iblock)
893 source_depth = float(sz / km)
894 conf_s = copy.deepcopy(self.qseis_s_config)
895 conf_r = copy.deepcopy(self.qseis_r_config)
897 gf_directory = op.abspath(self.qseis_baseconf.gf_directory)
899 fk_path = op.join(gf_directory, 'green_%.3fkm.fk' % source_depth)
900 info_path = op.join(gf_directory, 'green_%.3fkm.info' % source_depth)
902 conf_s.fk_path = fk_path
903 conf_s.info_path = info_path
905 conf_r.fk_path = fk_path
906 conf_r.info_path = info_path
908 if self.step == 0 and os.path.isfile(fk_path):
909 logger.info('Skipping step %i / %i, block %i / %i'
910 '(GF already exists)' %
911 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
912 return
914 logger.info(
915 'Starting step %i / %i, block %i / %i' %
916 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
918 dx = self.gf_config.distance_delta
919 conf_r.wavelet_duration = 0.001 * self.gf_config.sample_rate
921 if self.step == 0:
922 conf_s.source_depth = source_depth
923 runner = QSeisSRunner(tmp=self.tmp)
924 runner.run(conf_s)
926 else:
927 conf_r.receiver = QSeisRReceiver(lat=90 - firstx * cake.m2d,
928 lon=180.,
929 tstart=0.0,
930 distance=firstx)
931 conf_r.source = QSeis2dSource(lat=90 - 0.001 * dx * cake.m2d,
932 lon=0.0,
933 depth=source_depth)
935 runner = QSeisRRunner(tmp=self.tmp)
937 mmt1 = (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)),
938 {'r': (0, 1), 't': (3, 1), 'z': (5, 1)})
939 mmt2 = (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)),
940 {'r': (1, 1), 't': (4, 1), 'z': (6, 1)})
941 mmt3 = (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)),
942 {'r': (2, 1), 'z': (7, 1)})
943 mmt4 = (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)),
944 {'r': (8, 1), 'z': (9, 1)})
946 gfmapping = [mmt1, mmt2, mmt3, mmt4]
948 for mt, gfmap in gfmapping:
949 if mt:
950 m = mt.m()
951 f = float
952 conf_r.source_mech = QSeisRSourceMechMT(
953 mnn=f(m[0, 0]), mee=f(m[1, 1]), mdd=f(m[2, 2]),
954 mne=f(m[0, 1]), mnd=f(m[0, 2]), med=f(m[1, 2]))
955 else:
956 conf_r.source_mech = None
958 if conf_r.source_mech is not None:
959 runner.run(conf_r)
961 rawtraces = runner.get_traces()
963 interrupted = []
965 def signal_handler(signum, frame):
966 interrupted.append(True)
968 original = signal.signal(signal.SIGINT, signal_handler)
969 self.store.lock()
970 duplicate_inserts = 0
971 try:
972 for itr, tr in enumerate(rawtraces):
973 if tr.channel not in gfmap:
974 continue
976 x = tr.meta['distance']
977 if x > firstx + (nx - 1) * dx:
978 continue
980 ig, factor = gfmap[tr.channel]
982 if len(self.store.config.ns) == 2:
983 args = (sz, x, ig)
984 else:
985 args = (rz, sz, x, ig)
987 if self.qseis_baseconf.cut:
988 tmin = self.store.t(
989 self.qseis_baseconf.cut[0], args[:-1])
990 tmax = self.store.t(
991 self.qseis_baseconf.cut[1], args[:-1])
993 if None in (tmin, tmax):
994 continue
996 tr.chop(tmin, tmax)
998 tmin = tr.tmin
999 tmax = tr.tmax
1001 if self.qseis_baseconf.fade:
1002 ta, tb, tc, td = [
1003 self.store.t(v, args[:-1])
1004 for v in self.qseis_baseconf.fade]
1006 if None in (ta, tb, tc, td):
1007 continue
1009 if not (ta <= tb and tb <= tc and tc <= td):
1010 raise QSeis2dError(
1011 'invalid fade configuration')
1013 t = tr.get_xdata()
1014 fin = num.interp(t, [ta, tb], [0., 1.])
1015 fout = num.interp(t, [tc, td], [1., 0.])
1016 anti_fin = 1. - fin
1017 anti_fout = 1. - fout
1019 y = tr.ydata
1021 sum_anti_fin = num.sum(anti_fin)
1022 sum_anti_fout = num.sum(anti_fout)
1024 if sum_anti_fin != 0.0:
1025 yin = num.sum(anti_fin * y) / sum_anti_fin
1026 else:
1027 yin = 0.0
1029 if sum_anti_fout != 0.0:
1030 yout = num.sum(anti_fout * y) / sum_anti_fout
1031 else:
1032 yout = 0.0
1034 y2 = anti_fin * yin + \
1035 fin * fout * y + \
1036 anti_fout * yout
1038 if self.qseis_baseconf.relevel_with_fade_in:
1039 y2 -= yin
1041 tr.set_ydata(y2)
1043 gf_tr = gf.store.GFTrace.from_trace(tr)
1044 gf_tr.data *= factor
1046 try:
1047 self.store.put(args, gf_tr)
1048 except gf.store.DuplicateInsert:
1049 duplicate_inserts += 1
1051 finally:
1052 if duplicate_inserts:
1053 logger.warning('%i insertions skipped (duplicates)' %
1054 duplicate_inserts)
1056 self.store.unlock()
1057 signal.signal(signal.SIGINT, original)
1059 if interrupted:
1060 raise KeyboardInterrupt()
1062 logger.info(
1063 'Done with step %i / %i, block %i / %i' %
1064 (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
1067def init(store_dir, variant):
1068 if variant is None:
1069 variant = '2014'
1071 if variant not in ('2014'):
1072 raise gf.store.StoreError('unsupported variant: %s' % variant)
1074 modelling_code_id = 'qseis2d.%s' % variant
1076 qseis2d = QSeis2dConfig()
1078 qseis2d.time_region = (
1079 gf.meta.Timing('begin-50'),
1080 gf.meta.Timing('end+100'))
1082 qseis2d.cut = (
1083 gf.meta.Timing('begin-50'),
1084 gf.meta.Timing('end+100'))
1086 qseis2d.qseis_s_config.sw_flat_earth_transform = 1
1088 store_id = os.path.basename(os.path.realpath(store_dir))
1090 config = gf.meta.ConfigTypeA(
1092 id=store_id,
1093 ncomponents=10,
1094 sample_rate=0.2,
1095 receiver_depth=0 * km,
1096 source_depth_min=10 * km,
1097 source_depth_max=20 * km,
1098 source_depth_delta=10 * km,
1099 distance_min=100 * km,
1100 distance_max=1000 * km,
1101 distance_delta=10 * km,
1102 earthmodel_1d=cake.load_model().extract(depth_max='cmb'),
1103 modelling_code_id=modelling_code_id,
1104 tabulated_phases=[
1105 gf.meta.TPDef(
1106 id='begin',
1107 definition='p,P,p\\,P\\,Pv_(cmb)p'),
1108 gf.meta.TPDef(
1109 id='end',
1110 definition='2.5'),
1111 gf.meta.TPDef(
1112 id='P',
1113 definition='!P'),
1114 gf.meta.TPDef(
1115 id='S',
1116 definition='!S'),
1117 gf.meta.TPDef(
1118 id='p',
1119 definition='!p'),
1120 gf.meta.TPDef(
1121 id='s',
1122 definition='!s')])
1124 config.validate()
1125 return gf.store.Store.create_editables(
1126 store_dir, config=config, extra={'qseis2d': qseis2d})
1129def build(store_dir, force=False, nworkers=None, continue_=False, step=None,
1130 iblock=None):
1132 return QSeis2dGFBuilder.build(
1133 store_dir, force=force, nworkers=nworkers, continue_=continue_,
1134 step=step, iblock=iblock)