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

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 

17try: 

18 from scipy.integrate import cumulative_trapezoid 

19except ImportError: 

20 from scipy.integrate import cumtrapz as cumulative_trapezoid 

21 

22from pyrocko.moment_tensor import MomentTensor, symmat6 

23from pyrocko.guts import Float, Int, Tuple, List, Bool, Object, String 

24from pyrocko import trace, util, cake, gf 

25 

26km = 1e3 

27 

28guts_prefix = 'pf' 

29 

30Timing = gf.meta.Timing 

31 

32 

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

34 

35# how to call the programs 

36program_bins = { 

37 'qseis2d.qseisS2014': 'fomosto_qseisS2014', 

38 'qseis2d.qseisR2014': 'fomosto_qseisR2014', 

39} 

40 

41qseis2d_components = { 

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

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

44} 

45 

46# defaults 

47default_gf_directory = 'qseis2d_green' 

48default_fk_basefilename = 'green' 

49default_source_depth = 10.0 

50default_time_region = (Timing('-10'), Timing('+890')) 

51default_slowness_window = (0.0, 0.0, 0.2, 0.25) 

52 

53 

54def have_backend(): 

55 for cmd in [[exe] for exe in program_bins.values()]: 

56 try: 

57 p = Popen(cmd, stdout=PIPE, stderr=PIPE, stdin=PIPE) 

58 (stdout, stderr) = p.communicate() 

59 

60 except OSError: 

61 return False 

62 

63 return True 

64 

65 

66def nextpow2(i): 

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

68 

69 

70def str_float_vals(vals): 

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

72 

73 

74def str_int_vals(vals): 

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

76 

77 

78def str_str_vals(vals): 

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

80 

81 

82def scl(cs): 

83 if not cs: 

84 return '\n#' 

85 

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

87 

88 

89def cake_model_to_config(mod): 

90 k = 1000. 

91 srows = [] 

92 for i, row in enumerate(mod.to_scanlines()): 

93 depth, vp, vs, rho, qp, qs = row 

94 row = [depth / k, vp / k, vs / k, rho / k, qp, qs] 

95 srows.append('%i %s' % (i + 1, str_float_vals(row))) 

96 

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

98 

99 

100class QSeis2dSource(Object): 

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

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

103 depth = Float.T(default=default_source_depth) 

104 

105 def string_for_config(self): 

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

107 

108 

109class QSeisRSourceMech(Object): 

110 pass 

111 

112 

113class QSeisRSourceMechMT(QSeisRSourceMech): 

114 mnn = Float.T(default=1.0) 

115 mee = Float.T(default=1.0) 

116 mdd = Float.T(default=1.0) 

117 mne = Float.T(default=0.0) 

118 mnd = Float.T(default=0.0) 

119 med = Float.T(default=0.0) 

120 

121 def string_for_config(self): 

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

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

124 

125 

126class QSeisSPropagationFilter(Object): 

127 min_depth = Float.T(default=0.0) 

128 max_depth = Float.T(default=0.0) 

129 filtered_phase = Int.T(default=0) 

130 

131 def string_for_config(self): 

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

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

134 

135 

136class QSeisRBandpassFilter(Object): 

137 order = Int.T(default=-1) 

138 lower_cutoff = Float.T(default=0.1) # [Hz] 

139 upper_cutoff = Float.T(default=10.0) 

140 

141 def string_for_config(self): 

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

143 

144 

145class QSeisSConfig(Object): 

146 

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

148 

149 calc_slowness_window = Int.T(default=1) 

150 slowness_window = Tuple.T(4, optional=True) 

151 wavenumber_sampling = Float.T(default=2.5) 

152 aliasing_suppression_factor = Float.T(default=0.01) 

153 

154 filter_shallow_paths = Int.T(default=0) 

155 filter_shallow_paths_depth = Float.T(default=0.0) 

156 propagation_filters = List.T(QSeisSPropagationFilter.T()) 

157 

158 sw_flat_earth_transform = Int.T(default=0) 

159 

160 gradient_resolution_vp = Float.T(default=0.0) 

161 gradient_resolution_vs = Float.T(default=0.0) 

162 gradient_resolution_density = Float.T(default=0.0) 

163 

164 def items(self): 

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

166 

167 

168class QSeisSConfigFull(QSeisSConfig): 

169 

170 time_window = Float.T(default=900.0) 

171 

172 source_depth = Float.T(default=10.0) 

173 

174 receiver_basement_depth = Float.T(default=35.0) # [km] 

175 receiver_min_distance = Float.T(default=1000.0) # [km] 

176 receiver_max_distance = Float.T(default=10000.0) # [km] 

177 nsamples = Int.T(default=256) 

178 

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

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

181 

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

183 

184 @staticmethod 

185 def example(): 

186 conf = QSeisSConfigFull() 

187 conf.source_depth = 15. 

188 conf.receiver_basement_depth = 35. 

189 conf.receiver_max_distance = 2000. 

190 conf.earthmodel_1d = cake.load_model().extract(depth_max='cmb') 

191 conf.sw_flat_earth_transform = 1 

192 return conf 

193 

194 def string_for_config(self): 

195 

196 def aggregate(xx): 

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

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

199 

200 assert self.earthmodel_1d is not None 

201 assert self.slowness_window is not None or self.calc_slowness_window, \ 

202 "'slowness window' undefined and 'calc_slowness_window'=0" 

203 

204 d = self.__dict__.copy() 

205 

206 model_str, nlines = cake_model_to_config(self.earthmodel_1d) 

207 d['n_model_lines'] = nlines 

208 d['model_lines'] = model_str 

209 

210 if not self.slowness_window: 

211 d['str_slowness_window'] = str_float_vals(default_slowness_window) 

212 else: 

213 d['str_slowness_window'] = str_float_vals(self.slowness_window) 

214 

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

216 aggregate(self.propagation_filters) 

217 

218 template = '''# autogenerated QSEISS input by qseis2d.py 

219# This is the input file of FORTRAN77 program "QseisS" for calculation of 

220# f-k spectra of upgoing seismic waves at the reveiver-site basement. 

221# 

222# by 

223# Rongjiang Wang <wang@gfz-potsdam.de> 

224# GeoForschungsZentrum Potsdam 

225# Telegrafenberg, D-14473 Potsdam, Germany 

226# 

227# Modified from qseis2006, Potsdam, Dec. 2014 

228# 

229# SOURCE DEPTH 

230# ============ 

231# 1. source depth [km] 

232#------------------------------------------------------------------------------ 

233 %(source_depth)e |dble; 

234#------------------------------------------------------------------------------ 

235# 

236# RECEIVER-SITE PARAMETERS 

237# ======================== 

238# 1. receiver-site basement depth [km] 

239# 2. max. epicental distance [km] 

240#------------------------------------------------------------------------------ 

241 %(receiver_basement_depth)e |dble; 

242 %(receiver_max_distance)e |dble; 

243#------------------------------------------------------------------------------ 

244# TIME SAMPLING PARAMETERS 

245# ======================== 

246# 1. length of time window [sec] 

247# 2. number of time samples (<= 2*nfmax in qsglobal.h) 

248#------------------------------------------------------------------------------ 

249 %(time_window)e |dble: t_window; 

250 %(nsamples)i |int: no_t_samples; 

251#------------------------------------------------------------------------------ 

252# SLOWNESS WINDOW PARAMETERS 

253# ========================== 

254# 1. the low and high slowness cut-offs [s/km] with tapering: 0 < slw1 < slw2 

255# defining the bounds of cosine taper at the lower end, and 0 < slw3 < slw4 

256# defining the bounds of cosine taper at the higher end. 

257# 2. slowness sampling rate (1.0 = sampling with the Nyquist slowness, 2.0 = 

258# sampling with twice higher than the Nyquist slowness, and so on: the 

259# larger this parameter, the smaller the space-domain aliasing effect, but 

260# also the more computation effort); 

261# 3. the factor for suppressing time domain aliasing (> 0 and <= 1). 

262#------------------------------------------------------------------------------ 

263 %(str_slowness_window)s |dble: slw(1-4); 

264 %(wavenumber_sampling)e |dble: sample_rate; 

265 %(aliasing_suppression_factor)e |dble: supp_factor; 

266#------------------------------------------------------------------------------ 

267# OPTIONS FOR PARTIAL SOLUTIONS 

268# ============================= 

269# 1. switch for filtering waves with a shallow penetration depth (concerning 

270# their whole trace from source to receiver), penetration depth limit [km] 

271# (Note: if this option is selected, waves whose travel path never exceeds 

272# the given depth limit will be filtered ("seismic nuting"). the condition 

273# for selecting this filter is that the given shallow path depth limit 

274# should be larger than both source and receiver depth.) 

275# 2. number of depth ranges where the following selected up/down-sp2oing P or 

276# SV waves should be filtered 

277# 3. the 1. depth range: upper and lower depth [km], switch for filtering P 

278# or SV wave in this depth range: 

279# switch no: 1 2 3 4 other 

280# filtered phase: P(up) P(down) SV(up) SV(down) Error 

281# 4. the 2. ... 

282# (Note: the partial solution options are useful tools to increase the 

283# numerical resolution of desired wave phases. especially when the desired 

284# phases are much smaller than the undesired phases, these options should 

285# be selected and carefully combined.) 

286#------------------------------------------------------------------------------ 

287 %(filter_shallow_paths)i %(filter_shallow_paths_depth)e 

288 %(n_depth_ranges)i %(str_depth_ranges)s |int, dble, dble, int; 

289#------------------------------------------------------------------------------ 

290# OUTPUT FILE FOR F-K SPECTRA of GREEN'S FUNCTIONS 

291# ================================================ 

292# 1. info-file of Green's functions (ascii including f-k sampling parameters) 

293# 2. file name of Green's functions (binary files including explosion, strike 

294# -slip, dip-slip and clvd sources) 

295#------------------------------------------------------------------------------ 

296 '%(info_path)s' 

297 '%(fk_path)s' 

298#------------------------------------------------------------------------------ 

299# GLOBAL MODEL PARAMETERS 

300# ======================= 

301# 1. switch for flat-earth-transform 

302# 2. gradient resolution [percent] of vp, vs, and rho (density), if <= 0, then 

303# default values (depending on wave length at cut-off frequency) will be 

304# used 

305#------------------------------------------------------------------------------ 

306 %(sw_flat_earth_transform)i |int: sw_flat_earth_transform; 

307 %(gradient_resolution_vp)e %(gradient_resolution_vs)e %(gradient_resolution_density)e 

308#------------------------------------------------------------------------------ 

309# SOURCE-SITE LAYERED EARTH MODEL 

310# =============================== 

311# 1. number of data lines of the layered model 

312#------------------------------------------------------------------------------ 

313 %(n_model_lines)i |int: no_model_lines; 

314#------------------------------------------------------------------------------ 

315# no depth[km] vp[km/s] vs[km/s] rho[g/cm^3] qp qs 

316#------------------------------------------------------------------------------ 

317%(model_lines)s ; 

318#-----------------------END OF INPUT PARAMETERS-------------------------------- 

319 

320Glossary 

321 

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

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

324 

325SUPPRESSION OF TIME-DOMAIN ALIASING: The suppression of the time domain aliasing 

326is achieved by using the complex frequency technique. The suppression factor 

327should be a value between 0 and 1. If this factor is set to 0.1, for example, the 

328aliasing phase at the reduced time begin is suppressed to 10 percent. 

329 

330MODEL PARAMETER GRADIENT RESOLUTION: Layers with a constant gradient will be 

331discretized with a number of homogeneous sublayers. The gradient resolutions are 

332then used to determine the maximum allowed thickness of the sublayers. If the 

333resolutions of Vp, Vs and Rho (density) require different thicknesses, the 

334smallest is first chosen. If this is even smaller than 1 percent of the 

335characteristic wavelength, then the latter is taken finally for the sublayer 

336thickness. 

337''' # noqa 

338 return (template % d).encode('ascii') 

339 

340 

341class QSeisRReceiver(Object): 

342 lat = Float.T(default=10.0) 

343 lon = Float.T(default=0.0) 

344 depth = Float.T(default=0.0) 

345 tstart = Float.T(default=0.0) 

346 distance = Float.T(default=0.0) 

347 

348 def string_for_config(self): 

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

350 

351 

352class QSeisRConfig(Object): 

353 

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

355 

356 receiver_filter = QSeisRBandpassFilter.T(optional=True) 

357 

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

359 wavelet_type = Int.T(default=1) 

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

361 

362 def items(self): 

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

364 

365 

366class QSeisRConfigFull(QSeisRConfig): 

367 

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

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

370 

371 source_mech = QSeisRSourceMech.T( 

372 optional=True, 

373 default=QSeisRSourceMechMT.D()) 

374 

375 time_reduction = Float.T(default=0.0) 

376 

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

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

379 

380 output_format = Int.T(default=1) # 1/2 components in [Z, R, T]/[E, N, U] 

381 output_filename = String.T(default='seis.dat') 

382 

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

384 

385 @staticmethod 

386 def example(): 

387 conf = QSeisRConfigFull() 

388 conf.source = QSeis2dSource(lat=-80.5, lon=90.1) 

389 conf.receiver_location = QSeisRReceiver(lat=13.4, lon=240.5, depth=0.0) 

390 conf.time_reduction = 10.0 

391 conf.earthmodel_receiver_1d = cake.load_model().extract( 

392 depth_max='moho') 

393 return conf 

394 

395 @property 

396 def components(self): 

397 return qseis2d_components[self.output_format] 

398 

399 def get_output_filename(self, rundir): 

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

401 

402 def string_for_config(self): 

403 

404 def aggregate(xx): 

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

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

407 

408 assert self.earthmodel_receiver_1d is not None 

409 

410 d = self.__dict__.copy() 

411 

412# Actually not existing anymore in code 

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

414 

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

416 

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

418 

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

420 

421 model_str, nlines = cake_model_to_config(self.earthmodel_receiver_1d) 

422 

423 d['n_model_receiver_lines'] = nlines 

424 d['model_receiver_lines'] = model_str 

425 

426 if self.wavelet_type == 0: # user wavelet 

427 d['str_w_samples'] = '\n' \ 

428 + '%i\n' % len(self.user_wavelet_samples) \ 

429 + str_float_vals(self.user_wavelet_samples) 

430 else: 

431 d['str_w_samples'] = '' 

432 

433 if self.receiver_filter: 

434 d['str_receiver_filter'] = self.receiver_filter.string_for_config() 

435 else: 

436 d['str_receiver_filter'] = '-1 0.0 200.' 

437 

438 if self.source_mech: 

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

440 else: 

441 d['str_source'] = '0' 

442 

443 template = '''# autogenerated QSEISR input by qseis2d.py 

444# This is the input file of FORTRAN77 program "QseisR" for calculation of 

445# synthetic seismograms using the given f-k spectra of incident seismic 

446# waves at the receiver-site basement, where the f-k spectra should be 

447# prepared by the "QseisS" code. 

448# 

449# by 

450# Rongjiang Wang <wang@gfz-potsdam.de> 

451# GeoForschungsZentrum Potsdam 

452# Telegrafenberg, D-14473 Potsdam, Germany 

453# 

454# Modified from qseis2006, Potsdam, Dec. 2014 

455# 

456#------------------------------------------------------------------------------ 

457# SOURCE PARAMETERS 

458# ================= 

459# 1. epicenter (lat[deg], lon[deg]) 

460# 2. moment tensor in N*m: Mxx, Myy, Mzz, Mxy, Myz, Mzx 

461# Note: x is northward, y is eastward and z is downard 

462# conversion from CMT system: 

463# Mxx = Mtt, Myy= Mpp, Mzz = Mrr, Mxy = -Mtp, Myz = -Mrp, Mzx = Mrt 

464# 3. file of f-k spectra of incident waves 

465# 4. source duration [s], and selection of source time functions, i.e., 

466# source wavelet (0 = user's own wavelet; 1 = default wavelet: normalized 

467# square half-sinusoid) 

468# 5. if user's own wavelet is selected, then number of the wavelet time samples 

469# (<= 1024), followed by 

470# 6. the equidistant wavelet time series (no comment lines between the time 

471# series) 

472#------------------------------------------------------------------------------ 

473 %(str_source_location)s |dble(2); 

474 %(str_source)s |dble(6); 

475 '%(fk_path)s' |char; 

476 %(wavelet_duration)e %(wavelet_type)i %(str_w_samples)s |dble, int, dbls; 

477#------------------------------------------------------------------------------ 

478# RECEIVER PARAMETERS 

479# =================== 

480# 1. station location (lat[deg], lon[deg], depth[km]) 

481# Note: the epicentral distance should not exceed the max. distance used 

482# for generating the f-k spectra 

483# 2. time reduction[s] 

484# 3. order of bandpass filter(if <= 0, then no filter will be used), lower 

485# and upper cutoff frequences[Hz] 

486# 4. selection of output format (1 = Down/Radial/Azimuthal, 2 = East/North/Up) 

487# 5. output file of velocity seismograms 

488# 6. number of datalines representing the layered receiver-site structure, and 

489# selection of surface condition (0 = free surface, 1 = without free 

490# surface, i.e., upper halfspace is replaced by medium of the 1. layer) 

491# 7. ... layered structure model 

492#------------------------------------------------------------------------------ 

493 %(str_receiver)s |dble(3); 

494 %(time_reduction)e |dble; 

495 %(str_receiver_filter)s |int, dble(2); 

496 %(output_format)i |int; 

497 %(str_output_filename)s |char; 

498#------------------------------------------------------------------------------ 

499 %(n_model_receiver_lines)i |int: no_model_lines; 

500#------------------------------------------------------------------------------ 

501# MULTILAYERED MODEL PARAMETERS (shallow receiver-site structure) 

502# =============================================================== 

503# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs 

504#------------------------------------------------------------------------------ 

505%(model_receiver_lines)s 

506#-----------------------END OF INPUT PARAMETERS-------------------------------- 

507# 

508#-Requirements to use QSEIS2d: ------------------------------------------------ 

509# (1) Teleseismic body waves with penetration depth much larger than the 

510# receiver-site basement depth 

511# (2) The last layer parameters of the receiver-site structure should be 

512# identical with that of the source-site model at the depth which is 

513# defined as the common basement depth 

514# (3) The cutoff frequency should be high enough for separating different 

515# wave types. 

516''' # noqa 

517 return (template % d).encode('ascii') 

518 

519 

520class QSeis2dConfig(Object): 

521 ''' 

522 Combined config object for QSeisS and QSeisR. 

523 

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

525 the backend-independant Pyrocko GF Store config. 

526 ''' 

527 

528 qseis_s_config = QSeisSConfig.T(default=QSeisSConfig.D()) 

529 qseis_r_config = QSeisRConfig.T(default=QSeisRConfig.D()) 

530 qseis2d_version = String.T(default='2014') 

531 

532 time_region = Tuple.T(2, Timing.T(), default=default_time_region) 

533 cut = Tuple.T(2, Timing.T(), optional=True) 

534 fade = Tuple.T(4, Timing.T(), optional=True) 

535 relevel_with_fade_in = Bool.T(default=False) 

536 

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

538 

539 

540class QSeis2dError(gf.store.StoreError): 

541 pass 

542 

543 

544class Interrupted(gf.store.StoreError): 

545 def __str__(self): 

546 return 'Interrupted.' 

547 

548 

549class QSeisSRunner(object): 

550 ''' 

551 Takes QSeis2dConfigFull or QSeisSConfigFull objects, runs the program. 

552 ''' 

553 def __init__(self, tmp, keep_tmp=False): 

554 self.tempdir = mkdtemp(prefix='qseisSrun-', dir=tmp) 

555 self.keep_tmp = keep_tmp 

556 self.config = None 

557 

558 def run(self, config): 

559 

560 if isinstance(config, QSeis2dConfig): 

561 config = QSeisSConfig(**config.qseis_s_config) 

562 

563 self.config = config 

564 

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

566 

567 with open(input_fn, 'wb') as f: 

568 input_str = config.string_for_config() 

569 logger.debug('===== begin qseisS input =====\n' 

570 '%s===== end qseisS input =====' % input_str.decode()) 

571 f.write(input_str) 

572 

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

574 

575 old_wd = os.getcwd() 

576 os.chdir(self.tempdir) 

577 

578 interrupted = [] 

579 

580 def signal_handler(signum, frame): 

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

582 interrupted.append(True) 

583 

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

585 try: 

586 try: 

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

588 

589 except OSError: 

590 os.chdir(old_wd) 

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

592 

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

594 

595 finally: 

596 signal.signal(signal.SIGINT, original) 

597 

598 if interrupted: 

599 raise KeyboardInterrupt() 

600 

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

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

603 

604 errmess = [] 

605 if proc.returncode != 0: 

606 errmess.append( 

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

608 

609 if error_str: 

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

611 

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

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

614 

615 if errmess: 

616 self.keep_tmp = True 

617 

618 os.chdir(old_wd) 

619 raise QSeis2dError(''' 

620===== begin qseisS input ===== 

621%s===== end qseisS input ===== 

622===== begin qseisS output ===== 

623%s===== end qseisS output ===== 

624===== begin qseisS error ===== 

625%s===== end qseisS error ===== 

626%s 

627qseisS has been invoked as "%s" 

628in the directory %s'''.lstrip() % ( 

629 input_str.decode(), 

630 output_str.decode(), 

631 error_str.decode(), 

632 '\n'.join(errmess), program, 

633 self.tempdir)) 

634 

635 self.qseiss_output = output_str 

636 self.qseiss_error = error_str 

637 

638 os.chdir(old_wd) 

639 

640 def __del__(self): 

641 if self.tempdir: 

642 if not self.keep_tmp: 

643 shutil.rmtree(self.tempdir) 

644 self.tempdir = None 

645 else: 

646 logger.warning( 

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

648 

649 

650class QSeisRRunner(object): 

651 ''' 

652 Takes QSeis2dConfig or QSeisRConfigFull objects, runs the program and 

653 reads the output. 

654 ''' 

655 def __init__(self, tmp=None, keep_tmp=False): 

656 self.tempdir = mkdtemp(prefix='qseisRrun-', dir=tmp) 

657 self.keep_tmp = keep_tmp 

658 self.config = None 

659 

660 def run(self, config): 

661 if isinstance(config, QSeis2dConfig): 

662 config = QSeisRConfigFull(**config.qseis_r_config) 

663 

664 self.config = config 

665 

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

667 

668 with open(input_fn, 'wb') as f: 

669 input_str = config.string_for_config() 

670 old_wd = os.getcwd() 

671 os.chdir(self.tempdir) 

672 logger.debug('===== begin qseisR input =====\n' 

673 '%s===== end qseisR input =====' % input_str.decode()) 

674 f.write(input_str) 

675 

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

677 

678 interrupted = [] 

679 

680 def signal_handler(signum, frame): 

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

682 interrupted.append(True) 

683 

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

685 try: 

686 try: 

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

688 

689 except OSError: 

690 os.chdir(old_wd) 

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

692 

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

694 

695 finally: 

696 signal.signal(signal.SIGINT, original) 

697 

698 if interrupted: 

699 raise KeyboardInterrupt() 

700 

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

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

703 

704 errmess = [] 

705 if proc.returncode != 0: 

706 errmess.append( 

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

708 

709 if error_str: 

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

711 

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

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

714 

715 if errmess: 

716 self.keep_tmp = True 

717 

718 os.chdir(old_wd) 

719 raise QSeis2dError(''' 

720===== begin qseisR input ===== 

721%s===== end qseisR input ===== 

722===== begin qseisR output ===== 

723%s===== end qseisR output ===== 

724===== begin qseisR error ===== 

725%s===== end qseisR error ===== 

726%s 

727qseisR has been invoked as "%s" 

728in the directory %s'''.lstrip() % ( 

729 input_str.decode(), 

730 output_str.decode(), 

731 error_str.decode(), 

732 '\n'.join(errmess), program, 

733 self.tempdir)) 

734 

735 self.qseisr_output = output_str 

736 self.qseisr_error = error_str 

737 

738 os.chdir(old_wd) 

739 

740 def get_traces(self): 

741 

742 fn = self.config.get_output_filename(self.tempdir) 

743 data = num.loadtxt(fn, skiprows=1, dtype=float) 

744 nsamples, ntraces = data.shape 

745 deltat = (data[-1, 0] - data[0, 0]) / (nsamples - 1) 

746 toffset = data[0, 0] 

747 

748 tred = self.config.time_reduction 

749 rec = self.config.receiver 

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

751 

752 traces = [] 

753 for itrace, comp in enumerate(self.config.components): 

754 # qseis2d gives velocity-integrate to displacement 

755 # integration removes one sample, add it again in front 

756 displ = cumulative_trapezoid(num.concatenate( 

757 (num.zeros(1), data[:, itrace + 1])), dx=deltat) 

758 

759 tr = trace.Trace( 

760 '', '%04i' % itrace, '', comp, 

761 tmin=tmin, deltat=deltat, ydata=displ, 

762 meta=dict(distance=rec.distance, 

763 azimuth=0.0)) 

764 

765 traces.append(tr) 

766 

767 return traces 

768 

769 def __del__(self): 

770 if self.tempdir: 

771 if not self.keep_tmp: 

772 shutil.rmtree(self.tempdir) 

773 self.tempdir = None 

774 else: 

775 logger.warning( 

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

777 

778 

779class QSeis2dGFBuilder(gf.builder.Builder): 

780 nsteps = 2 

781 

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

783 force=False): 

784 

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

786 

787 storeconf = self.store.config 

788 

789 if step == 0: 

790 block_size = (1, 1, storeconf.ndistances) 

791 else: 

792 if block_size is None: 

793 block_size = (1, 1, 1) # QSeisR does only allow one receiver 

794 

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

796 block_size = block_size[1:] 

797 

798 gf.builder.Builder.__init__( 

799 self, storeconf, step, block_size=block_size) 

800 

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

802 

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

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

805 

806 conf_s.earthmodel_1d = storeconf.earthmodel_1d 

807 if storeconf.earthmodel_receiver_1d is not None: 

808 conf_r.earthmodel_receiver_1d = \ 

809 storeconf.earthmodel_receiver_1d 

810 

811 else: 

812 conf_r.earthmodel_receiver_1d = \ 

813 storeconf.earthmodel_1d.extract( 

814 depth_max='moho') 

815 # depth_max=conf_s.receiver_basement_depth*km) 

816 

817 deltat = 1.0 / self.gf_config.sample_rate 

818 

819 if 'time_window_min' not in shared: 

820 d = self.store.make_timing_params( 

821 baseconf.time_region[0], baseconf.time_region[1], 

822 force=force) 

823 

824 shared['time_window_min'] = float( 

825 num.ceil(d['tlenmax'] / self.gf_config.sample_rate) * 

826 self.gf_config.sample_rate) 

827 shared['time_reduction'] = d['tmin_vred'] 

828 

829 time_window_min = shared['time_window_min'] 

830 

831 conf_s.nsamples = nextpow2(int(round(time_window_min / deltat)) + 1) 

832 conf_s.time_window = (conf_s.nsamples - 1) * deltat 

833 conf_r.time_reduction = shared['time_reduction'] 

834 

835 if step == 0: 

836 if 'slowness_window' not in shared: 

837 if conf_s.calc_slowness_window: 

838 phases = [ 

839 storeconf.tabulated_phases[i].phases 

840 for i in range(len( 

841 storeconf.tabulated_phases))] 

842 

843 all_phases = [] 

844 map(all_phases.extend, phases) 

845 

846 mean_source_depth = num.mean(( 

847 storeconf.source_depth_min, 

848 storeconf.source_depth_max)) 

849 

850 arrivals = conf_s.earthmodel_1d.arrivals( 

851 phases=all_phases, 

852 distances=num.linspace( 

853 conf_s.receiver_min_distance, 

854 conf_s.receiver_max_distance, 

855 100) * cake.m2d, 

856 zstart=mean_source_depth) 

857 

858 ps = num.array( 

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

860 

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

862 

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

864 0., 

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

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

867 

868 else: 

869 shared['slowness_window'] = conf_s.slowness_window