Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/fomosto/qseis2d.py: 79%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6import numpy as num 

7import logging 

8import os 

9import shutil 

10import math 

11import copy 

12import signal 

13 

14from tempfile import mkdtemp 

15from subprocess import Popen, PIPE 

16import os.path as op 

17from scipy.integrate import cumtrapz 

18 

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 

22 

23km = 1e3 

24 

25guts_prefix = 'pf' 

26 

27Timing = gf.meta.Timing 

28 

29 

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

31 

32# how to call the programs 

33program_bins = { 

34 'qseis2d.qseisS2014': 'fomosto_qseisS2014', 

35 'qseis2d.qseisR2014': 'fomosto_qseisR2014', 

36} 

37 

38qseis2d_components = { 

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

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

41} 

42 

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) 

49 

50 

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

56 

57 except OSError: 

58 return False 

59 

60 return True 

61 

62 

63def nextpow2(i): 

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

65 

66 

67def str_float_vals(vals): 

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

69 

70 

71def str_int_vals(vals): 

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

73 

74 

75def str_str_vals(vals): 

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

77 

78 

79def scl(cs): 

80 if not cs: 

81 return '\n#' 

82 

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

84 

85 

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

93 

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

95 

96 

97class QSeis2dSource(Object): 

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

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

100 depth = Float.T(default=default_source_depth) 

101 

102 def string_for_config(self): 

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

104 

105 

106class QSeisRSourceMech(Object): 

107 pass 

108 

109 

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) 

117 

118 def string_for_config(self): 

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

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

121 

122 

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) 

127 

128 def string_for_config(self): 

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

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

131 

132 

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) 

137 

138 def string_for_config(self): 

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

140 

141 

142class QSeisSConfig(Object): 

143 

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

145 

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) 

150 

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

154 

155 sw_flat_earth_transform = Int.T(default=0) 

156 

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) 

160 

161 def items(self): 

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

163 

164 

165class QSeisSConfigFull(QSeisSConfig): 

166 

167 time_window = Float.T(default=900.0) 

168 

169 source_depth = Float.T(default=10.0) 

170 

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) 

175 

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

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

178 

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

180 

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 

190 

191 def string_for_config(self): 

192 

193 def aggregate(xx): 

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

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

196 

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" 

200 

201 d = self.__dict__.copy() 

202 

203 model_str, nlines = cake_model_to_config(self.earthmodel_1d) 

204 d['n_model_lines'] = nlines 

205 d['model_lines'] = model_str 

206 

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) 

211 

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

213 aggregate(self.propagation_filters) 

214 

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

316 

317Glossary 

318 

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

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

321 

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. 

326 

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

336 

337 

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) 

344 

345 def string_for_config(self): 

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

347 

348 

349class QSeisRConfig(Object): 

350 

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

352 

353 receiver_filter = QSeisRBandpassFilter.T(optional=True) 

354 

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

356 wavelet_type = Int.T(default=1) 

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

358 

359 def items(self): 

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

361 

362 

363class QSeisRConfigFull(QSeisRConfig): 

364 

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

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

367 

368 source_mech = QSeisRSourceMech.T( 

369 optional=True, 

370 default=QSeisRSourceMechMT.D()) 

371 

372 time_reduction = Float.T(default=0.0) 

373 

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

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

376 

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

379 

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

381 

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 

391 

392 @property 

393 def components(self): 

394 return qseis2d_components[self.output_format] 

395 

396 def get_output_filename(self, rundir): 

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

398 

399 def string_for_config(self): 

400 

401 def aggregate(xx): 

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

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

404 

405 assert self.earthmodel_receiver_1d is not None 

406 

407 d = self.__dict__.copy() 

408 

409# Actually not existing anymore in code 

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

411 

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

413 

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

415 

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

417 

418 model_str, nlines = cake_model_to_config(self.earthmodel_receiver_1d) 

419 

420 d['n_model_receiver_lines'] = nlines 

421 d['model_receiver_lines'] = model_str 

422 

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

429 

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

434 

435 if self.source_mech: 

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

437 else: 

438 d['str_source'] = '0' 

439 

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

515 

516 

517class QSeis2dConfig(Object): 

518 ''' 

519 Combined config object for QSeisS and QSeisR. 

520 

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

522 the backend-independant Pyrocko GF Store config. 

523 ''' 

524 

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

528 

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) 

533 

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

535 

536 

537class QSeis2dError(gf.store.StoreError): 

538 pass 

539 

540 

541class Interrupted(gf.store.StoreError): 

542 def __str__(self): 

543 return 'Interrupted.' 

544 

545 

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 

554 

555 def run(self, config): 

556 

557 if isinstance(config, QSeis2dConfig): 

558 config = QSeisSConfig(**config.qseis_s_config) 

559 

560 self.config = config 

561 

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

563 

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) 

569 

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

571 

572 old_wd = os.getcwd() 

573 os.chdir(self.tempdir) 

574 

575 interrupted = [] 

576 

577 def signal_handler(signum, frame): 

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

579 interrupted.append(True) 

580 

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

582 try: 

583 try: 

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

585 

586 except OSError: 

587 os.chdir(old_wd) 

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

589 

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

591 

592 finally: 

593 signal.signal(signal.SIGINT, original) 

594 

595 if interrupted: 

596 raise KeyboardInterrupt() 

597 

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

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

600 

601 errmess = [] 

602 if proc.returncode != 0: 

603 errmess.append( 

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

605 

606 if error_str: 

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

608 

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

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

611 

612 if errmess: 

613 self.keep_tmp = True 

614 

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

631 

632 self.qseiss_output = output_str 

633 self.qseiss_error = error_str 

634 

635 os.chdir(old_wd) 

636 

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) 

645 

646 

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 

656 

657 def run(self, config): 

658 if isinstance(config, QSeis2dConfig): 

659 config = QSeisRConfigFull(**config.qseis_r_config) 

660 

661 self.config = config 

662 

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

664 

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) 

672 

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

674 

675 interrupted = [] 

676 

677 def signal_handler(signum, frame): 

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

679 interrupted.append(True) 

680 

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

682 try: 

683 try: 

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

685 

686 except OSError: 

687 os.chdir(old_wd) 

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

689 

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

691 

692 finally: 

693 signal.signal(signal.SIGINT, original) 

694 

695 if interrupted: 

696 raise KeyboardInterrupt() 

697 

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

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

700 

701 errmess = [] 

702 if proc.returncode != 0: 

703 errmess.append( 

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

705 

706 if error_str: 

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

708 

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

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

711 

712 if errmess: 

713 self.keep_tmp = True 

714 

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

731 

732 self.qseisr_output = output_str 

733 self.qseisr_error = error_str 

734 

735 os.chdir(old_wd) 

736 

737 def get_traces(self): 

738 

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] 

744 

745 tred = self.config.time_reduction 

746 rec = self.config.receiver 

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

748 

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) 

755 

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

761 

762 traces.append(tr) 

763 

764 return traces 

765 

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) 

774 

775 

776class QSeis2dGFBuilder(gf.builder.Builder): 

777 nsteps = 2 

778 

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

780 force=False): 

781 

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

783 

784 storeconf = self.store.config 

785 

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 

791 

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

793 block_size = block_size[1:] 

794 

795 gf.builder.Builder.__init__( 

796 self, storeconf, step, block_size=block_size) 

797 

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

799 

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

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

802 

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 

807 

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) 

813 

814 deltat = 1.0 / self.gf_config.sample_rate 

815 

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) 

820 

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

825 

826 time_window_min = shared['time_window_min'] 

827 

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

831 

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

839 

840 all_phases = [] 

841 map(all_phases.extend, phases) 

842 

843 mean_source_depth = num.mean(( 

844 storeconf.source_depth_min, 

845 storeconf.source_depth_max)) 

846 

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) 

854 

855 ps = num.array( 

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

857 

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

859 

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

861 0., 

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

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

864 

865 else: 

866 shared['slowness_window'] = conf_s.slowness_window 

867 

868 conf_s.slowness_window = shared['slowness_window'] 

869