1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5from __future__ import absolute_import, division 

6 

7import numpy as num 

8import logging 

9import os 

10import shutil 

11import math 

12import copy 

13import signal 

14 

15from tempfile import mkdtemp 

16from subprocess import Popen, PIPE 

17import os.path as op 

18from scipy.integrate import cumtrapz 

19 

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 

23 

24km = 1e3 

25 

26guts_prefix = 'pf' 

27 

28Timing = gf.meta.Timing 

29 

30 

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

32 

33# how to call the programs 

34program_bins = { 

35 'qseis2d.qseisS2014': 'fomosto_qseisS2014', 

36 'qseis2d.qseisR2014': 'fomosto_qseisR2014', 

37} 

38 

39qseis2d_components = { 

40 1: 'z r t'.split(), 

41 2: 'e n u'.split(), 

42} 

43 

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) 

50 

51 

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() 

57 

58 except OSError: 

59 return False 

60 

61 return True 

62 

63 

64def nextpow2(i): 

65 return 2 ** int(math.ceil(math.log(i) / math.log(2.))) 

66 

67 

68def str_float_vals(vals): 

69 return ' '.join('%e' % val for val in vals) 

70 

71 

72def str_int_vals(vals): 

73 return ' '.join('%i' % val for val in vals) 

74 

75 

76def str_str_vals(vals): 

77 return ' '.join("'%s'" % val for val in vals) 

78 

79 

80def scl(cs): 

81 if not cs: 

82 return '\n#' 

83 

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

85 

86 

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))) 

94 

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

96 

97 

98class QSeis2dSource(Object): 

99 lat = Float.T(default=10.) 

100 lon = Float.T(default=0.) 

101 depth = Float.T(default=default_source_depth) 

102 

103 def string_for_config(self): 

104 return '%(lat)e %(lon)15e ' % self.__dict__ 

105 

106 

107class QSeisRSourceMech(Object): 

108 pass 

109 

110 

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) 

118 

119 def string_for_config(self): 

120 return '%(mnn)e %(mee)15e %(mdd)15e ' \ 

121 '%(mne)15e %(med)15e %(mnd)15e ' % self.__dict__ 

122 

123 

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) 

128 

129 def string_for_config(self): 

130 return '%(min_depth)15e %(max_depth)15e ' \ 

131 '%(filtered_phase)i' % self.__dict__ 

132 

133 

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) 

138 

139 def string_for_config(self): 

140 return '%(order) %(lower_cutoff)5e %(upper_cutoff)5e' % self.__dict__ 

141 

142 

143class QSeisSConfig(Object): 

144 

145 qseiss_version = String.T(default='2014') 

146 

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) 

151 

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()) 

155 

156 sw_flat_earth_transform = Int.T(default=0) 

157 

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) 

161 

162 def items(self): 

163 return dict(self.T.inamevals(self)) 

164 

165 

166class QSeisSConfigFull(QSeisSConfig): 

167 

168 time_window = Float.T(default=900.0) 

169 

170 source_depth = Float.T(default=10.0) 

171 

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) 

176 

177 info_path = String.T(default='green.info') 

178 fk_path = String.T(default='green.fk') 

179 

180 earthmodel_1d = gf.meta.Earthmodel1D.T(optional=True) 

181 

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 

191 

192 def string_for_config(self): 

193 

194 def aggregate(xx): 

195 return len(xx), '\n'.join( 

196 [''] + [x.string_for_config() for x in xx]) 

197 

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" 

201 

202 d = self.__dict__.copy() 

203 

204 model_str, nlines = cake_model_to_config(self.earthmodel_1d) 

205 d['n_model_lines'] = nlines 

206 d['model_lines'] = model_str 

207 

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) 

212 

213 d['n_depth_ranges'], d['str_depth_ranges'] = \ 

214 aggregate(self.propagation_filters) 

215 

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-------------------------------- 

317 

318Glossary 

319 

320SLOWNESS: The slowness is the inverse of apparent wave velocity = sin(i)/v with 

321i = incident angle and v = true wave velocity. 

322 

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. 

327 

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') 

337 

338 

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) 

345 

346 def string_for_config(self): 

347 return '%(lat)e %(lon)15e %(depth)15e' % self.__dict__ 

348 

349 

350class QSeisRConfig(Object): 

351 

352 qseisr_version = String.T(default='2014') 

353 

354 receiver_filter = QSeisRBandpassFilter.T(optional=True) 

355 

356 wavelet_duration = Float.T(default=0.001) # [s] 

357 wavelet_type = Int.T(default=1) 

358 user_wavelet_samples = List.T(Float.T()) 

359 

360 def items(self): 

361 return dict(self.T.inamevals(self)) 

362 

363 

364class QSeisRConfigFull(QSeisRConfig): 

365 

366 source = QSeis2dSource(depth=default_source_depth) # [lat, lon(deg)] 

367 receiver = QSeisRReceiver() # [lat, lon(deg)] 

368 

369 source_mech = QSeisRSourceMech.T( 

370 optional=True, 

371 default=QSeisRSourceMechMT.D()) 

372 

373 time_reduction = Float.T(default=0.0) 

374 

375 info_path = String.T(default='green.info') 

376 fk_path = String.T(default='green.fk') 

377 

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') 

380 

381 earthmodel_receiver_1d = gf.meta.Earthmodel1D.T(optional=True) 

382 

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 

392 

393 @property 

394 def components(self): 

395 return qseis2d_components[self.output_format] 

396 

397 def get_output_filename(self, rundir): 

398 return op.join(rundir, self.output_filename) 

399 

400 def string_for_config(self): 

401 

402 def aggregate(xx): 

403 return len(xx), '\n'.join( 

404 [''] + [x.string_for_config() for x in xx]) 

405 

406 assert self.earthmodel_receiver_1d is not None 

407 

408 d = self.__dict__.copy() 

409 

410# Actually not existing anymore in code 

411# d['sw_surface'] = 0 # 0-free-surface, 1-no fre surface 

412 

413 d['str_source_location'] = self.source.string_for_config() 

414 

415 d['str_receiver'] = self.receiver.string_for_config() 

416 

417 d['str_output_filename'] = "'%s'" % self.output_filename 

418 

419 model_str, nlines = cake_model_to_config(self.earthmodel_receiver_1d) 

420 

421 d['n_model_receiver_lines'] = nlines 

422 d['model_receiver_lines'] = model_str 

423 

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'] = '' 

430 

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.' 

435 

436 if self.source_mech: 

437 d['str_source'] = '%s' % (self.source_mech.string_for_config()) 

438 else: 

439 d['str_source'] = '0' 

440 

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') 

516 

517 

518class QSeis2dConfig(Object): 

519 ''' 

520 Combined config object for QSeisS and QSeisR. 

521 

522 This config object should contain all settings which cannot be derived from 

523 the backend-independant Pyrocko GF Store config. 

524 ''' 

525 

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') 

529 

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) 

534 

535 gf_directory = String.T(default='qseis2d_green') 

536 

537 

538class QSeis2dError(gf.store.StoreError): 

539 pass 

540 

541 

542class Interrupted(gf.store.StoreError): 

543 def __str__(self): 

544 return 'Interrupted.' 

545 

546 

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 

555 

556 def run(self, config): 

557 

558 if isinstance(config, QSeis2dConfig): 

559 config = QSeisSConfig(**config.qseis_s_config) 

560 

561 self.config = config 

562 

563 input_fn = op.join(self.tempdir, 'input') 

564 

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) 

570 

571 program = program_bins['qseis2d.qseisS%s' % config.qseiss_version] 

572 

573 old_wd = os.getcwd() 

574 os.chdir(self.tempdir) 

575 

576 interrupted = [] 

577 

578 def signal_handler(signum, frame): 

579 os.kill(proc.pid, signal.SIGTERM) 

580 interrupted.append(True) 

581 

582 original = signal.signal(signal.SIGINT, signal_handler) 

583 try: 

584 try: 

585 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE) 

586 

587 except OSError: 

588 os.chdir(old_wd) 

589 raise QSeis2dError('could not start qseisS: "%s"' % program) 

590 

591 (output_str, error_str) = proc.communicate(b'input\n') 

592 

593 finally: 

594 signal.signal(signal.SIGINT, original) 

595 

596 if interrupted: 

597 raise KeyboardInterrupt() 

598 

599 logger.debug('===== begin qseisS output =====\n' 

600 '%s===== end qseisS output =====' % output_str.decode()) 

601 

602 errmess = [] 

603 if proc.returncode != 0: 

604 errmess.append( 

605 'qseisS had a non-zero exit state: %i' % proc.returncode) 

606 

607 if error_str: 

608 errmess.append('qseisS emitted something via stderr') 

609 

610 if output_str.lower().find(b'error') != -1: 

611 errmess.append("the string 'error' appeared in qseisS output") 

612 

613 if errmess: 

614 self.keep_tmp = True 

615 

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)) 

632 

633 self.qseiss_output = output_str 

634 self.qseiss_error = error_str 

635 

636 os.chdir(old_wd) 

637 

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.warning( 

645 'not removing temporary directory: %s' % self.tempdir) 

646 

647 

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 

657 

658 def run(self, config): 

659 if isinstance(config, QSeis2dConfig): 

660 config = QSeisRConfigFull(**config.qseis_r_config) 

661 

662 self.config = config 

663 

664 input_fn = op.join(self.tempdir, 'input') 

665 

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) 

673 

674 program = program_bins['qseis2d.qseisR%s' % config.qseisr_version] 

675 

676 interrupted = [] 

677 

678 def signal_handler(signum, frame): 

679 os.kill(proc.pid, signal.SIGTERM) 

680 interrupted.append(True) 

681 

682 original = signal.signal(signal.SIGINT, signal_handler) 

683 try: 

684 try: 

685 proc = Popen(program, stdin=PIPE, stdout=PIPE, stderr=PIPE) 

686 

687 except OSError: 

688 os.chdir(old_wd) 

689 raise QSeis2dError('could not start qseisR: "%s"' % program) 

690 

691 (output_str, error_str) = proc.communicate(b'input\n') 

692 

693 finally: 

694 signal.signal(signal.SIGINT, original) 

695 

696 if interrupted: 

697 raise KeyboardInterrupt() 

698 

699 logger.debug('===== begin qseisR output =====\n' 

700 '%s===== end qseisR output =====' % output_str.decode()) 

701 

702 errmess = [] 

703 if proc.returncode != 0: 

704 errmess.append( 

705 'qseisR had a non-zero exit state: %i' % proc.returncode) 

706 

707 if error_str: 

708 errmess.append('qseisR emitted something via stderr') 

709 

710 if output_str.lower().find(b'error') != -1: 

711 errmess.append("the string 'error' appeared in qseisR output") 

712 

713 if errmess: 

714 self.keep_tmp = True 

715 

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)) 

732 

733 self.qseisr_output = output_str 

734 self.qseisr_error = error_str 

735 

736 os.chdir(old_wd) 

737 

738 def get_traces(self): 

739 

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] 

745 

746 tred = self.config.time_reduction 

747 rec = self.config.receiver 

748 tmin = rec.tstart + toffset + deltat + tred 

749 

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) 

756 

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)) 

762 

763 traces.append(tr) 

764 

765 return traces 

766 

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.warning( 

774 'not removing temporary directory: %s' % self.tempdir) 

775 

776 

777class QSeis2dGFBuilder(gf.builder.Builder): 

778 nsteps = 2 

779 

780 def __init__(self, store_dir, step, shared, block_size=None, tmp=None, 

781 force=False): 

782 

783 self.store = gf.store.Store(store_dir, 'w') 

784 

785 storeconf = self.store.config 

786 

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 

792 

793 if len(storeconf.ns) == 2: 

794 block_size = block_size[1:] 

795 

796 gf.builder.Builder.__init__( 

797 self, storeconf, step, block_size=block_size) 

798 

799 baseconf = self.store.get_extra('qseis2d') 

800 

801 conf_s = QSeisSConfigFull(**baseconf.qseis_s_config.items()) 

802 conf_r = QSeisRConfigFull(**baseconf.qseis_r_config.items()) 

803 

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 

808 

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) 

814 

815 deltat = 1.0 / self.gf_config.sample_rate 

816 

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) 

821 

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'] 

826 

827 time_window_min = shared['time_window_min'] 

828 

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'] 

832 

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))] 

840 

841 all_phases = [] 

842 map(all_phases.extend, phases) 

843 

844 mean_source_depth = num.mean(( 

845 storeconf.source_depth_min, 

846 storeconf.source_depth_max)) 

847 

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) 

855 

856 ps = num.array( 

857 [arrivals[i].p for i in range(len(arrivals))]) 

858 

859 slownesses = ps / (cake.r2d * cake.d2m / km) 

860 

861 shared['slowness_window'] = (0., 

862 0., 

863 1.1 * float(slownesses.max()), 

864 1.3 * float(slownesses.max())) 

865 

866 else: 

867 shared['slowness_window'] = conf_s.slowness_window 

868 

869 conf_s.slowness_window = shared['slowness_window'] 

870 

871 self.qseis_s_config = conf_s 

872 self.qseis_r_config = conf_r 

873 self.qseis_baseconf = baseconf 

874 

875 self.tmp = tmp 

876 if self.tmp is not None: 

877 util.ensuredir(self.tmp) 

878 

879 util.ensuredir(baseconf.gf_directory) 

880 

881 def cleanup(self): 

882 self.store.close() 

883 

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) 

888 

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) 

893 

894 source_depth = float(sz / km) 

895 conf_s = copy.deepcopy(self.qseis_s_config) 

896 conf_r = copy.deepcopy(self.qseis_r_config) 

897 

898 gf_directory = op.abspath(self.qseis_baseconf.gf_directory) 

899 

900 fk_path = op.join(gf_directory, 'green_%.3fkm.fk' % source_depth) 

901 info_path = op.join(gf_directory, 'green_%.3fkm.info' % source_depth) 

902 

903 conf_s.fk_path = fk_path 

904 conf_s.info_path = info_path 

905 

906 conf_r.fk_path = fk_path 

907 conf_r.info_path = info_path 

908 

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 

914 

915 logger.info( 

916 'Starting step %i / %i, block %i / %i' % 

917 (self.step + 1, self.nsteps, iblock + 1, self.nblocks)) 

918 

919 dx = self.gf_config.distance_delta 

920 conf_r.wavelet_duration = 0.001 * self.gf_config.sample_rate 

921 

922 if self.step == 0: 

923 conf_s.source_depth = source_depth 

924 runner = QSeisSRunner(tmp=self.tmp) 

925 runner.run(conf_s) 

926 

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) 

935 

936 runner = QSeisRRunner(tmp=self.tmp) 

937 

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)}) 

946 

947 gfmapping = [mmt1, mmt2, mmt3, mmt4] 

948 

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 

958 

959 if conf_r.source_mech is not None: 

960 runner.run(conf_r) 

961 

962 rawtraces = runner.get_traces() 

963 

964 interrupted = [] 

965 

966 def signal_handler(signum, frame): 

967 interrupted.append(True) 

968 

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 

976 

977 x = tr.meta['distance'] 

978 if x > firstx + (nx - 1) * dx: 

979 continue 

980 

981 ig, factor = gfmap[tr.channel] 

982 

983 if len(self.store.config.ns) == 2: 

984 args = (sz, x, ig) 

985 else: 

986 args = (rz, sz, x, ig) 

987 

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]) 

993 

994 if None in (tmin, tmax): 

995 continue 

996 

997 tr.chop(tmin, tmax) 

998 

999 tmin = tr.tmin 

1000 tmax = tr.tmax 

1001 

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] 

1006 

1007 if None in (ta, tb, tc, td): 

1008 continue 

1009 

1010 if not (ta <= tb and tb <= tc and tc <= td): 

1011 raise QSeis2dError( 

1012 'invalid fade configuration') 

1013 

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 

1019 

1020 y = tr.ydata 

1021 

1022 sum_anti_fin = num.sum(anti_fin) 

1023 sum_anti_fout = num.sum(anti_fout) 

1024 

1025 if sum_anti_fin != 0.0: 

1026 yin = num.sum(anti_fin * y) / sum_anti_fin 

1027 else: 

1028 yin = 0.0 

1029 

1030 if sum_anti_fout != 0.0: 

1031 yout = num.sum(anti_fout * y) / sum_anti_fout 

1032 else: 

1033 yout = 0.0 

1034 

1035 y2 = anti_fin * yin + \ 

1036 fin * fout * y + \ 

1037 anti_fout * yout 

1038 

1039 if self.qseis_baseconf.relevel_with_fade_in: 

1040 y2 -= yin 

1041 

1042 tr.set_ydata(y2) 

1043 

1044 gf_tr = gf.store.GFTrace.from_trace(tr) 

1045 gf_tr.data *= factor 

1046 

1047 try: 

1048 self.store.put(args, gf_tr) 

1049 except gf.store.DuplicateInsert: 

1050 duplicate_inserts += 1 

1051 

1052 finally: 

1053 if duplicate_inserts: 

1054 logger.warning('%i insertions skipped (duplicates)' % 

1055 duplicate_inserts) 

1056 

1057 self.store.unlock() 

1058 signal.signal(signal.SIGINT, original) 

1059 

1060 if interrupted: 

1061 raise KeyboardInterrupt() 

1062 

1063 logger.info( 

1064 'Done with step %i / %i, block %i / %i' % 

1065 (self.step + 1, self.nsteps, iblock + 1, self.nblocks)) 

1066 

1067 

1068def init(store_dir, variant): 

1069 if variant is None: 

1070 variant = '2014' 

1071 

1072 if variant not in ('2014'): 

1073 raise gf.store.StoreError('unsupported variant: %s' % variant) 

1074 

1075 modelling_code_id = 'qseis2d.%s' % variant 

1076 

1077 qseis2d = QSeis2dConfig() 

1078 

1079 qseis2d.time_region = ( 

1080 gf.meta.Timing('begin-50'), 

1081 gf.meta.Timing('end+100')) 

1082 

1083 qseis2d.cut = ( 

1084 gf.meta.Timing('begin-50'), 

1085 gf.meta.Timing('end+100')) 

1086 

1087 qseis2d.qseis_s_config.sw_flat_earth_transform = 1 

1088 

1089 store_id = os.path.basename(os.path.realpath(store_dir)) 

1090 

1091 config = gf.meta.ConfigTypeA( 

1092 

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')]) 

1124 

1125 config.validate() 

1126 return gf.store.Store.create_editables( 

1127 store_dir, config=config, extra={'qseis2d': qseis2d}) 

1128 

1129 

1130def build(store_dir, force=False, nworkers=None, continue_=False, step=None, 

1131 iblock=None): 

1132 

1133 return QSeis2dGFBuilder.build( 

1134 store_dir, force=force, nworkers=nworkers, continue_=continue_, 

1135 step=step, iblock=iblock)