Source code for pyrocko.fomosto.qseis2d

# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------

import numpy as num
import logging
import os
import shutil
import math
import copy
import signal

from tempfile import mkdtemp
from subprocess import Popen, PIPE
import os.path as op
from scipy.integrate import cumtrapz

from pyrocko.moment_tensor import MomentTensor, symmat6
from pyrocko.guts import Float, Int, Tuple, List, Bool, Object, String
from pyrocko import trace, util, cake, gf

km = 1e3

guts_prefix = 'pf'

Timing = gf.meta.Timing


logger = logging.getLogger('pyrocko.fomosto.qseis2d')

# how to call the programs
program_bins = {
    'qseis2d.qseisS2014': 'fomosto_qseisS2014',
    'qseis2d.qseisR2014': 'fomosto_qseisR2014',
}

qseis2d_components = {
    1: 'z r t'.split(),
    2: 'e n u'.split(),
}

# defaults
default_gf_directory = 'qseis2d_green'
default_fk_basefilename = 'green'
default_source_depth = 10.0
default_time_region = (Timing('-10'), Timing('+890'))
default_slowness_window = (0.0, 0.0, 0.2, 0.25)


def have_backend():
    for cmd in [[exe] for exe in program_bins.values()]:
        try:
            p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE)
            (stdout, stderr) = p.communicate()

        except OSError:
            return False

    return True


def nextpow2(i):
    return 2 ** int(math.ceil(math.log(i) / math.log(2.)))


def str_float_vals(vals):
    return ' '.join('%e' % val for val in vals)


def str_int_vals(vals):
    return ' '.join('%i' % val for val in vals)


def str_str_vals(vals):
    return ' '.join("'%s'" % val for val in vals)


def scl(cs):
    if not cs:
        return '\n#'

    return '\n' + ' '.join('(%e,%e)' % (c.real, c.imag) for c in cs)


def cake_model_to_config(mod):
    k = 1000.
    srows = []
    for i, row in enumerate(mod.to_scanlines()):
        depth, vp, vs, rho, qp, qs = row
        row = [depth / k, vp / k, vs / k, rho / k, qp, qs]
        srows.append('%i %s' % (i + 1, str_float_vals(row)))

    return '\n'.join(srows), len(srows)


[docs]class QSeis2dSource(Object): lat = Float.T(default=10.) lon = Float.T(default=0.) depth = Float.T(default=default_source_depth) def string_for_config(self): return '%(lat)e %(lon)15e ' % self.__dict__
[docs]class QSeisRSourceMech(Object): pass
[docs]class QSeisRSourceMechMT(QSeisRSourceMech): mnn = Float.T(default=1.0) mee = Float.T(default=1.0) mdd = Float.T(default=1.0) mne = Float.T(default=0.0) mnd = Float.T(default=0.0) med = Float.T(default=0.0) def string_for_config(self): return '%(mnn)e %(mee)15e %(mdd)15e ' \ '%(mne)15e %(med)15e %(mnd)15e ' % self.__dict__
[docs]class QSeisSPropagationFilter(Object): min_depth = Float.T(default=0.0) max_depth = Float.T(default=0.0) filtered_phase = Int.T(default=0) def string_for_config(self): return '%(min_depth)15e %(max_depth)15e ' \ '%(filtered_phase)i' % self.__dict__
[docs]class QSeisRBandpassFilter(Object): order = Int.T(default=-1) lower_cutoff = Float.T(default=0.1) # [Hz] upper_cutoff = Float.T(default=10.0) def string_for_config(self): return '%(order)i %(lower_cutoff)5e %(upper_cutoff)5e' % self.__dict__
[docs]class QSeisSConfig(Object): qseiss_version = String.T(default='2014') calc_slowness_window = Int.T(default=1) slowness_window = Tuple.T(4, optional=True) wavenumber_sampling = Float.T(default=2.5) aliasing_suppression_factor = Float.T(default=0.01) filter_shallow_paths = Int.T(default=0) filter_shallow_paths_depth = Float.T(default=0.0) propagation_filters = List.T(QSeisSPropagationFilter.T()) sw_flat_earth_transform = Int.T(default=0) gradient_resolution_vp = Float.T(default=0.0) gradient_resolution_vs = Float.T(default=0.0) gradient_resolution_density = Float.T(default=0.0) def items(self): return dict(self.T.inamevals(self))
[docs]class QSeisSConfigFull(QSeisSConfig): time_window = Float.T(default=900.0) source_depth = Float.T(default=10.0) receiver_basement_depth = Float.T(default=35.0) # [km] receiver_min_distance = Float.T(default=1000.0) # [km] receiver_max_distance = Float.T(default=10000.0) # [km] nsamples = Int.T(default=256) info_path = String.T(default='green.info') fk_path = String.T(default='green.fk') earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True) @staticmethod def example(): conf = QSeisSConfigFull() conf.source_depth = 15. conf.receiver_basement_depth = 35. conf.receiver_max_distance = 2000. conf.earthmodel_1d = cake.load_model().extract(depth_max='cmb') conf.sw_flat_earth_transform = 1 return conf def string_for_config(self): def aggregate(xx): return len(xx), '\n'.join( [''] + [x.string_for_config() for x in xx]) assert self.earthmodel_1d is not None assert self.slowness_window is not None or self.calc_slowness_window,\ "'slowness window' undefined and 'calc_slowness_window'=0" d = self.__dict__.copy() model_str, nlines = cake_model_to_config(self.earthmodel_1d) d['n_model_lines'] = nlines d['model_lines'] = model_str if not self.slowness_window: d['str_slowness_window'] = str_float_vals(default_slowness_window) else: d['str_slowness_window'] = str_float_vals(self.slowness_window) d['n_depth_ranges'], d['str_depth_ranges'] = \ aggregate(self.propagation_filters) template = '''# autogenerated QSEISS input by qseis2d.py # This is the input file of FORTRAN77 program "QseisS" for calculation of # f-k spectra of upgoing seismic waves at the reveiver-site basement. # # by # Rongjiang Wang <wang@gfz-potsdam.de> # GeoForschungsZentrum Potsdam # Telegrafenberg, D-14473 Potsdam, Germany # # Modified from qseis2006, Potsdam, Dec. 2014 # # SOURCE DEPTH # ============ # 1. source depth [km] #------------------------------------------------------------------------------ %(source_depth)e |dble; #------------------------------------------------------------------------------ # # RECEIVER-SITE PARAMETERS # ======================== # 1. receiver-site basement depth [km] # 2. max. epicental distance [km] #------------------------------------------------------------------------------ %(receiver_basement_depth)e |dble; %(receiver_max_distance)e |dble; #------------------------------------------------------------------------------ # TIME SAMPLING PARAMETERS # ======================== # 1. length of time window [sec] # 2. number of time samples (<= 2*nfmax in qsglobal.h) #------------------------------------------------------------------------------ %(time_window)e |dble: t_window; %(nsamples)i |int: no_t_samples; #------------------------------------------------------------------------------ # SLOWNESS WINDOW PARAMETERS # ========================== # 1. the low and high slowness cut-offs [s/km] with tapering: 0 < slw1 < slw2 # defining the bounds of cosine taper at the lower end, and 0 < slw3 < slw4 # defining the bounds of cosine taper at the higher end. # 2. slowness sampling rate (1.0 = sampling with the Nyquist slowness, 2.0 = # sampling with twice higher than the Nyquist slowness, and so on: the # larger this parameter, the smaller the space-domain aliasing effect, but # also the more computation effort); # 3. the factor for suppressing time domain aliasing (> 0 and <= 1). #------------------------------------------------------------------------------ %(str_slowness_window)s |dble: slw(1-4); %(wavenumber_sampling)e |dble: sample_rate; %(aliasing_suppression_factor)e |dble: supp_factor; #------------------------------------------------------------------------------ # OPTIONS FOR PARTIAL SOLUTIONS # ============================= # 1. switch for filtering waves with a shallow penetration depth (concerning # their whole trace from source to receiver), penetration depth limit [km] # (Note: if this option is selected, waves whose travel path never exceeds # the given depth limit will be filtered ("seismic nuting"). the condition # for selecting this filter is that the given shallow path depth limit # should be larger than both source and receiver depth.) # 2. number of depth ranges where the following selected up/down-sp2oing P or # SV waves should be filtered # 3. the 1. depth range: upper and lower depth [km], switch for filtering P # or SV wave in this depth range: # switch no: 1 2 3 4 other # filtered phase: P(up) P(down) SV(up) SV(down) Error # 4. the 2. ... # (Note: the partial solution options are useful tools to increase the # numerical resolution of desired wave phases. especially when the desired # phases are much smaller than the undesired phases, these options should # be selected and carefully combined.) #------------------------------------------------------------------------------ %(filter_shallow_paths)i %(filter_shallow_paths_depth)e %(n_depth_ranges)i %(str_depth_ranges)s |int, dble, dble, int; #------------------------------------------------------------------------------ # OUTPUT FILE FOR F-K SPECTRA of GREEN'S FUNCTIONS # ================================================ # 1. info-file of Green's functions (ascii including f-k sampling parameters) # 2. file name of Green's functions (binary files including explosion, strike # -slip, dip-slip and clvd sources) #------------------------------------------------------------------------------ '%(info_path)s' '%(fk_path)s' #------------------------------------------------------------------------------ # GLOBAL MODEL PARAMETERS # ======================= # 1. switch for flat-earth-transform # 2. gradient resolution [percent] of vp, vs, and rho (density), if <= 0, then # default values (depending on wave length at cut-off frequency) will be # used #------------------------------------------------------------------------------ %(sw_flat_earth_transform)i |int: sw_flat_earth_transform; %(gradient_resolution_vp)e %(gradient_resolution_vs)e %(gradient_resolution_density)e #------------------------------------------------------------------------------ # SOURCE-SITE LAYERED EARTH MODEL # =============================== # 1. number of data lines of the layered model #------------------------------------------------------------------------------ %(n_model_lines)i |int: no_model_lines; #------------------------------------------------------------------------------ # no depth[km] vp[km/s] vs[km/s] rho[g/cm^3] qp qs #------------------------------------------------------------------------------ %(model_lines)s ; #-----------------------END OF INPUT PARAMETERS-------------------------------- Glossary SLOWNESS: The slowness is the inverse of apparent wave velocity = sin(i)/v with i = incident angle and v = true wave velocity. SUPPRESSION OF TIME-DOMAIN ALIASING: The suppression of the time domain aliasing is achieved by using the complex frequency technique. The suppression factor should be a value between 0 and 1. If this factor is set to 0.1, for example, the aliasing phase at the reduced time begin is suppressed to 10 percent. MODEL PARAMETER GRADIENT RESOLUTION: Layers with a constant gradient will be discretized with a number of homogeneous sublayers. The gradient resolutions are then used to determine the maximum allowed thickness of the sublayers. If the resolutions of Vp, Vs and Rho (density) require different thicknesses, the smallest is first chosen. If this is even smaller than 1 percent of the characteristic wavelength, then the latter is taken finally for the sublayer thickness. ''' # noqa return (template % d).encode('ascii')
[docs]class QSeisRReceiver(Object): lat = Float.T(default=10.0) lon = Float.T(default=0.0) depth = Float.T(default=0.0) tstart = Float.T(default=0.0) distance = Float.T(default=0.0) def string_for_config(self): return '%(lat)e %(lon)15e %(depth)15e' % self.__dict__
[docs]class QSeisRConfig(Object): qseisr_version = String.T(default='2014') receiver_filter = QSeisRBandpassFilter.T(optional=True) wavelet_duration = Float.T(default=0.001) # [s] wavelet_type = Int.T(default=1) user_wavelet_samples = List.T(Float.T()) def items(self): return dict(self.T.inamevals(self))
[docs]class QSeisRConfigFull(QSeisRConfig): source = QSeis2dSource(depth=default_source_depth) # [lat, lon(deg)] receiver = QSeisRReceiver() # [lat, lon(deg)] source_mech = QSeisRSourceMech.T( optional=True, default=QSeisRSourceMechMT.D()) time_reduction = Float.T(default=0.0) info_path = String.T(default='green.info') fk_path = String.T(default='green.fk') output_format = Int.T(default=1) # 1/2 components in [Z, R, T]/[E, N, U] output_filename = String.T(default='seis.dat') earthmodel_receiver_1d = gf.meta.Earthmodel1D.T(optional=True) @staticmethod def example(): conf = QSeisRConfigFull() conf.source = QSeis2dSource(lat=-80.5, lon=90.1) conf.receiver_location = QSeisRReceiver(lat=13.4, lon=240.5, depth=0.0) conf.time_reduction = 10.0 conf.earthmodel_receiver_1d = cake.load_model().extract( depth_max='moho') return conf @property def components(self): return qseis2d_components[self.output_format] def get_output_filename(self, rundir): return op.join(rundir, self.output_filename) def string_for_config(self): def aggregate(xx): return len(xx), '\n'.join( [''] + [x.string_for_config() for x in xx]) assert self.earthmodel_receiver_1d is not None d = self.__dict__.copy() # Actually not existing anymore in code # d['sw_surface'] = 0 # 0-free-surface, 1-no fre surface d['str_source_location'] = self.source.string_for_config() d['str_receiver'] = self.receiver.string_for_config() d['str_output_filename'] = "'%s'" % self.output_filename model_str, nlines = cake_model_to_config(self.earthmodel_receiver_1d) d['n_model_receiver_lines'] = nlines d['model_receiver_lines'] = model_str if self.wavelet_type == 0: # user wavelet d['str_w_samples'] = '\n' \ + '%i\n' % len(self.user_wavelet_samples) \ + str_float_vals(self.user_wavelet_samples) else: d['str_w_samples'] = '' if self.receiver_filter: d['str_receiver_filter'] = self.receiver_filter.string_for_config() else: d['str_receiver_filter'] = '-1 0.0 200.' if self.source_mech: d['str_source'] = '%s' % (self.source_mech.string_for_config()) else: d['str_source'] = '0' template = '''# autogenerated QSEISR input by qseis2d.py # This is the input file of FORTRAN77 program "QseisR" for calculation of # synthetic seismograms using the given f-k spectra of incident seismic # waves at the receiver-site basement, where the f-k spectra should be # prepared by the "QseisS" code. # # by # Rongjiang Wang <wang@gfz-potsdam.de> # GeoForschungsZentrum Potsdam # Telegrafenberg, D-14473 Potsdam, Germany # # Modified from qseis2006, Potsdam, Dec. 2014 # #------------------------------------------------------------------------------ # SOURCE PARAMETERS # ================= # 1. epicenter (lat[deg], lon[deg]) # 2. moment tensor in N*m: Mxx, Myy, Mzz, Mxy, Myz, Mzx # Note: x is northward, y is eastward and z is downard # conversion from CMT system: # Mxx = Mtt, Myy= Mpp, Mzz = Mrr, Mxy = -Mtp, Myz = -Mrp, Mzx = Mrt # 3. file of f-k spectra of incident waves # 4. source duration [s], and selection of source time functions, i.e., # source wavelet (0 = user's own wavelet; 1 = default wavelet: normalized # square half-sinusoid) # 5. if user's own wavelet is selected, then number of the wavelet time samples # (<= 1024), followed by # 6. the equidistant wavelet time series (no comment lines between the time # series) #------------------------------------------------------------------------------ %(str_source_location)s |dble(2); %(str_source)s |dble(6); '%(fk_path)s' |char; %(wavelet_duration)e %(wavelet_type)i %(str_w_samples)s |dble, int, dbls; #------------------------------------------------------------------------------ # RECEIVER PARAMETERS # =================== # 1. station location (lat[deg], lon[deg], depth[km]) # Note: the epicentral distance should not exceed the max. distance used # for generating the f-k spectra # 2. time reduction[s] # 3. order of bandpass filter(if <= 0, then no filter will be used), lower # and upper cutoff frequences[Hz] # 4. selection of output format (1 = Down/Radial/Azimuthal, 2 = East/North/Up) # 5. output file of velocity seismograms # 6. number of datalines representing the layered receiver-site structure, and # selection of surface condition (0 = free surface, 1 = without free # surface, i.e., upper halfspace is replaced by medium of the 1. layer) # 7. ... layered structure model #------------------------------------------------------------------------------ %(str_receiver)s |dble(3); %(time_reduction)e |dble; %(str_receiver_filter)s |int, dble(2); %(output_format)i |int; %(str_output_filename)s |char; #------------------------------------------------------------------------------ %(n_model_receiver_lines)i |int: no_model_lines; #------------------------------------------------------------------------------ # MULTILAYERED MODEL PARAMETERS (shallow receiver-site structure) # =============================================================== # no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs #------------------------------------------------------------------------------ %(model_receiver_lines)s #-----------------------END OF INPUT PARAMETERS-------------------------------- # #-Requirements to use QSEIS2d: ------------------------------------------------ # (1) Teleseismic body waves with penetration depth much larger than the # receiver-site basement depth # (2) The last layer parameters of the receiver-site structure should be # identical with that of the source-site model at the depth which is # defined as the common basement depth # (3) The cutoff frequency should be high enough for separating different # wave types. ''' # noqa return (template % d).encode('ascii')
[docs]class QSeis2dConfig(Object): ''' Combined config object for QSeisS and QSeisR. This config object should contain all settings which cannot be derived from the backend-independant Pyrocko GF Store config. ''' qseis_s_config = QSeisSConfig.T(default=QSeisSConfig.D()) qseis_r_config = QSeisRConfig.T(default=QSeisRConfig.D()) qseis2d_version = String.T(default='2014') time_region = Tuple.T(2, Timing.T(), default=default_time_region) cut = Tuple.T(2, Timing.T(), optional=True) fade = Tuple.T(4, Timing.T(), optional=True) relevel_with_fade_in = Bool.T(default=False) gf_directory = String.T(default='qseis2d_green')
class QSeis2dError(gf.store.StoreError): pass class Interrupted(gf.store.StoreError): def __str__(self): return 'Interrupted.'
[docs]class QSeisSRunner(object): ''' Takes QSeis2dConfigFull or QSeisSConfigFull objects, runs the program. ''' def __init__(self, tmp, keep_tmp=False): self.tempdir = mkdtemp(prefix='qseisSrun-', dir=tmp) self.keep_tmp = keep_tmp self.config = None def run(self, config): if isinstance(config, QSeis2dConfig): config = QSeisSConfig(**config.qseis_s_config) self.config = config input_fn = op.join(self.tempdir, 'input') with open(input_fn, 'wb') as f: input_str = config.string_for_config() logger.debug('===== begin qseisS input =====\n' '%s===== end qseisS input =====' % input_str.decode()) f.write(input_str) program = program_bins['qseis2d.qseisS%s' % config.qseiss_version] old_wd = os.getcwd() os.chdir(self.tempdir) interrupted = [] def signal_handler(signum, frame): os.kill(proc.pid, signal.SIGTERM) interrupted.append(True) original = signal.signal(signal.SIGINT, signal_handler) try: try: proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE) except OSError: os.chdir(old_wd) raise QSeis2dError('could not start qseisS: "%s"' % program) (output_str, error_str) = proc.communicate(b'input\n') finally: signal.signal(signal.SIGINT, original) if interrupted: raise KeyboardInterrupt() logger.debug('===== begin qseisS output =====\n' '%s===== end qseisS output =====' % output_str.decode()) errmess = [] if proc.returncode != 0: errmess.append( 'qseisS had a non-zero exit state: %i' % proc.returncode) if error_str: errmess.append('qseisS emitted something via stderr') if output_str.lower().find(b'error') != -1: errmess.append("the string 'error' appeared in qseisS output") if errmess: self.keep_tmp = True os.chdir(old_wd) raise QSeis2dError(''' ===== begin qseisS input ===== %s===== end qseisS input ===== ===== begin qseisS output ===== %s===== end qseisS output ===== ===== begin qseisS error ===== %s===== end qseisS error ===== %s qseisS has been invoked as "%s" in the directory %s'''.lstrip() % ( input_str.decode(), output_str.decode(), error_str.decode(), '\n'.join(errmess), program, self.tempdir)) self.qseiss_output = output_str self.qseiss_error = error_str os.chdir(old_wd) def __del__(self): if self.tempdir: if not self.keep_tmp: shutil.rmtree(self.tempdir) self.tempdir = None else: logger.warning( 'not removing temporary directory: %s' % self.tempdir)
[docs]class QSeisRRunner(object): ''' Takes QSeis2dConfig or QSeisRConfigFull objects, runs the program and reads the output. ''' def __init__(self, tmp=None, keep_tmp=False): self.tempdir = mkdtemp(prefix='qseisRrun-', dir=tmp) self.keep_tmp = keep_tmp self.config = None def run(self, config): if isinstance(config, QSeis2dConfig): config = QSeisRConfigFull(**config.qseis_r_config) self.config = config input_fn = op.join(self.tempdir, 'input') with open(input_fn, 'wb') as f: input_str = config.string_for_config() old_wd = os.getcwd() os.chdir(self.tempdir) logger.debug('===== begin qseisR input =====\n' '%s===== end qseisR input =====' % input_str.decode()) f.write(input_str) program = program_bins['qseis2d.qseisR%s' % config.qseisr_version] interrupted = [] def signal_handler(signum, frame): os.kill(proc.pid, signal.SIGTERM) interrupted.append(True) original = signal.signal(signal.SIGINT, signal_handler) try: try: proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE) except OSError: os.chdir(old_wd) raise QSeis2dError('could not start qseisR: "%s"' % program) (output_str, error_str) = proc.communicate(b'input\n') finally: signal.signal(signal.SIGINT, original) if interrupted: raise KeyboardInterrupt() logger.debug('===== begin qseisR output =====\n' '%s===== end qseisR output =====' % output_str.decode()) errmess = [] if proc.returncode != 0: errmess.append( 'qseisR had a non-zero exit state: %i' % proc.returncode) if error_str: errmess.append('qseisR emitted something via stderr') if output_str.lower().find(b'error') != -1: errmess.append("the string 'error' appeared in qseisR output") if errmess: self.keep_tmp = True os.chdir(old_wd) raise QSeis2dError(''' ===== begin qseisR input ===== %s===== end qseisR input ===== ===== begin qseisR output ===== %s===== end qseisR output ===== ===== begin qseisR error ===== %s===== end qseisR error ===== %s qseisR has been invoked as "%s" in the directory %s'''.lstrip() % ( input_str.decode(), output_str.decode(), error_str.decode(), '\n'.join(errmess), program, self.tempdir)) self.qseisr_output = output_str self.qseisr_error = error_str os.chdir(old_wd) def get_traces(self): fn = self.config.get_output_filename(self.tempdir) data = num.loadtxt(fn, skiprows=1, dtype=float) nsamples, ntraces = data.shape deltat = (data[-1, 0] - data[0, 0]) / (nsamples - 1) toffset = data[0, 0] tred = self.config.time_reduction rec = self.config.receiver tmin = rec.tstart + toffset + deltat + tred traces = [] for itrace, comp in enumerate(self.config.components): # qseis2d gives velocity-integrate to displacement # integration removes one sample, add it again in front displ = cumtrapz(num.concatenate( (num.zeros(1), data[:, itrace + 1])), dx=deltat) tr = trace.Trace( '', '%04i' % itrace, '', comp, tmin=tmin, deltat=deltat, ydata=displ, meta=dict(distance=rec.distance, azimuth=0.0)) traces.append(tr) return traces def __del__(self): if self.tempdir: if not self.keep_tmp: shutil.rmtree(self.tempdir) self.tempdir = None else: logger.warning( 'not removing temporary directory: %s' % self.tempdir)
class QSeis2dGFBuilder(gf.builder.Builder): nsteps = 2 def __init__(self, store_dir, step, shared, block_size=None, tmp=None, force=False): self.store = gf.store.Store(store_dir, 'w') storeconf = self.store.config if step == 0: block_size = (1, 1, storeconf.ndistances) else: if block_size is None: block_size = (1, 1, 1) # QSeisR does only allow one receiver if len(storeconf.ns) == 2: block_size = block_size[1:] gf.builder.Builder.__init__( self, storeconf, step, block_size=block_size) baseconf = self.store.get_extra('qseis2d') conf_s = QSeisSConfigFull(**baseconf.qseis_s_config.items()) conf_r = QSeisRConfigFull(**baseconf.qseis_r_config.items()) conf_s.earthmodel_1d = storeconf.earthmodel_1d if storeconf.earthmodel_receiver_1d is not None: conf_r.earthmodel_receiver_1d = \ storeconf.earthmodel_receiver_1d else: conf_r.earthmodel_receiver_1d = \ storeconf.earthmodel_1d.extract( depth_max='moho') # depth_max=conf_s.receiver_basement_depth*km) deltat = 1.0 / self.gf_config.sample_rate if 'time_window_min' not in shared: d = self.store.make_timing_params( baseconf.time_region[0], baseconf.time_region[1], force=force) shared['time_window_min'] = float( num.ceil(d['tlenmax'] / self.gf_config.sample_rate) * self.gf_config.sample_rate) shared['time_reduction'] = d['tmin_vred'] time_window_min = shared['time_window_min'] conf_s.nsamples = nextpow2(int(round(time_window_min / deltat)) + 1) conf_s.time_window = (conf_s.nsamples - 1) * deltat conf_r.time_reduction = shared['time_reduction'] if step == 0: if 'slowness_window' not in shared: if conf_s.calc_slowness_window: phases = [ storeconf.tabulated_phases[i].phases for i in range(len( storeconf.tabulated_phases))] all_phases = [] map(all_phases.extend, phases) mean_source_depth = num.mean(( storeconf.source_depth_min, storeconf.source_depth_max)) arrivals = conf_s.earthmodel_1d.arrivals( phases=all_phases, distances=num.linspace( conf_s.receiver_min_distance, conf_s.receiver_max_distance, 100) * cake.m2d, zstart=mean_source_depth) ps = num.array( [arrivals[i].p for i in range(len(arrivals))]) slownesses = ps / (cake.r2d * cake.d2m / km) shared['slowness_window'] = (0., 0., 1.1 * float(slownesses.max()), 1.3 * float(slownesses.max())) else: shared['slowness_window'] = conf_s.slowness_window conf_s.slowness_window = shared['slowness_window'] self.qseis_s_config = conf_s self.qseis_r_config = conf_r self.qseis_baseconf = baseconf self.tmp = tmp if self.tmp is not None: util.ensuredir(self.tmp) util.ensuredir(baseconf.gf_directory) def cleanup(self): self.store.close() def work_block(self, iblock): if len(self.store.config.ns) == 2: (sz, firstx), (sz, lastx), (ns, nx) = \ self.get_block_extents(iblock) rz = self.store.config.receiver_depth else: (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \ self.get_block_extents(iblock) source_depth = float(sz / km) conf_s = copy.deepcopy(self.qseis_s_config) conf_r = copy.deepcopy(self.qseis_r_config) gf_directory = op.abspath(self.qseis_baseconf.gf_directory) fk_path = op.join(gf_directory, 'green_%.3fkm.fk' % source_depth) info_path = op.join(gf_directory, 'green_%.3fkm.info' % source_depth) conf_s.fk_path = fk_path conf_s.info_path = info_path conf_r.fk_path = fk_path conf_r.info_path = info_path if self.step == 0 and os.path.isfile(fk_path): logger.info('Skipping step %i / %i, block %i / %i' '(GF already exists)' % (self.step + 1, self.nsteps, iblock + 1, self.nblocks)) return logger.info( 'Starting step %i / %i, block %i / %i' % (self.step + 1, self.nsteps, iblock + 1, self.nblocks)) dx = self.gf_config.distance_delta conf_r.wavelet_duration = 0.001 * self.gf_config.sample_rate if self.step == 0: conf_s.source_depth = source_depth runner = QSeisSRunner(tmp=self.tmp) runner.run(conf_s) else: conf_r.receiver = QSeisRReceiver(lat=90 - firstx * cake.m2d, lon=180., tstart=0.0, distance=firstx) conf_r.source = QSeis2dSource(lat=90 - 0.001 * dx * cake.m2d, lon=0.0, depth=source_depth) runner = QSeisRRunner(tmp=self.tmp) mmt1 = (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)), {'r': (0, 1), 't': (3, 1), 'z': (5, 1)}) mmt2 = (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)), {'r': (1, 1), 't': (4, 1), 'z': (6, 1)}) mmt3 = (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)), {'r': (2, 1), 'z': (7, 1)}) mmt4 = (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)), {'r': (8, 1), 'z': (9, 1)}) gfmapping = [mmt1, mmt2, mmt3, mmt4] for mt, gfmap in gfmapping: if mt: m = mt.m() f = float conf_r.source_mech = QSeisRSourceMechMT( mnn=f(m[0, 0]), mee=f(m[1, 1]), mdd=f(m[2, 2]), mne=f(m[0, 1]), mnd=f(m[0, 2]), med=f(m[1, 2])) else: conf_r.source_mech = None if conf_r.source_mech is not None: runner.run(conf_r) rawtraces = runner.get_traces() interrupted = [] def signal_handler(signum, frame): interrupted.append(True) original = signal.signal(signal.SIGINT, signal_handler) self.store.lock() duplicate_inserts = 0 try: for itr, tr in enumerate(rawtraces): if tr.channel not in gfmap: continue x = tr.meta['distance'] if x > firstx + (nx - 1) * dx: continue ig, factor = gfmap[tr.channel] if len(self.store.config.ns) == 2: args = (sz, x, ig) else: args = (rz, sz, x, ig) if self.qseis_baseconf.cut: tmin = self.store.t( self.qseis_baseconf.cut[0], args[:-1]) tmax = self.store.t( self.qseis_baseconf.cut[1], args[:-1]) if None in (tmin, tmax): continue tr.chop(tmin, tmax) tmin = tr.tmin tmax = tr.tmax if self.qseis_baseconf.fade: ta, tb, tc, td = [ self.store.t(v, args[:-1]) for v in self.qseis_baseconf.fade] if None in (ta, tb, tc, td): continue if not (ta <= tb and tb <= tc and tc <= td): raise QSeis2dError( 'invalid fade configuration') t = tr.get_xdata() fin = num.interp(t, [ta, tb], [0., 1.]) fout = num.interp(t, [tc, td], [1., 0.]) anti_fin = 1. - fin anti_fout = 1. - fout y = tr.ydata sum_anti_fin = num.sum(anti_fin) sum_anti_fout = num.sum(anti_fout) if sum_anti_fin != 0.0: yin = num.sum(anti_fin * y) / sum_anti_fin else: yin = 0.0 if sum_anti_fout != 0.0: yout = num.sum(anti_fout * y) / sum_anti_fout else: yout = 0.0 y2 = anti_fin * yin + \ fin * fout * y + \ anti_fout * yout if self.qseis_baseconf.relevel_with_fade_in: y2 -= yin tr.set_ydata(y2) gf_tr = gf.store.GFTrace.from_trace(tr) gf_tr.data *= factor try: self.store.put(args, gf_tr) except gf.store.DuplicateInsert: duplicate_inserts += 1 finally: if duplicate_inserts: logger.warning('%i insertions skipped (duplicates)' % duplicate_inserts) self.store.unlock() signal.signal(signal.SIGINT, original) if interrupted: raise KeyboardInterrupt() logger.info( 'Done with step %i / %i, block %i / %i' % (self.step + 1, self.nsteps, iblock + 1, self.nblocks)) def init(store_dir, variant): if variant is None: variant = '2014' if variant not in ('2014'): raise gf.store.StoreError('unsupported variant: %s' % variant) modelling_code_id = 'qseis2d.%s' % variant qseis2d = QSeis2dConfig() qseis2d.time_region = ( gf.meta.Timing('begin-50'), gf.meta.Timing('end+100')) qseis2d.cut = ( gf.meta.Timing('begin-50'), gf.meta.Timing('end+100')) qseis2d.qseis_s_config.sw_flat_earth_transform = 1 store_id = os.path.basename(os.path.realpath(store_dir)) config = gf.meta.ConfigTypeA( id=store_id, ncomponents=10, sample_rate=0.2, receiver_depth=0 * km, source_depth_min=10 * km, source_depth_max=20 * km, source_depth_delta=10 * km, distance_min=100 * km, distance_max=1000 * km, distance_delta=10 * km, earthmodel_1d=cake.load_model().extract(depth_max='cmb'), modelling_code_id=modelling_code_id, tabulated_phases=[ gf.meta.TPDef( id='begin', definition='p,P,p\\,P\\,Pv_(cmb)p'), gf.meta.TPDef( id='end', definition='2.5'), gf.meta.TPDef( id='P', definition='!P'), gf.meta.TPDef( id='S', definition='!S'), gf.meta.TPDef( id='p', definition='!p'), gf.meta.TPDef( id='s', definition='!s')]) config.validate() return gf.store.Store.create_editables( store_dir, config=config, extra={'qseis2d': qseis2d}) def build(store_dir, force=False, nworkers=None, continue_=False, step=None, iblock=None): return QSeis2dGFBuilder.build( store_dir, force=force, nworkers=nworkers, continue_=continue_, step=step, iblock=iblock)