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

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 glob 

11import copy 

12import signal 

13import math 

14import time 

15 

16from tempfile import mkdtemp 

17from subprocess import Popen, PIPE 

18from os.path import join as pjoin 

19 

20from pyrocko import trace, util, cake, gf 

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

22from pyrocko.moment_tensor import MomentTensor, symmat6 

23 

24guts_prefix = 'pf' 

25 

26logger = logging.getLogger('pyrocko.fomosto.qssp') 

27 

28# how to call the programs 

29program_bins = { 

30 'qssp.2010beta': 'fomosto_qssp2010beta', 

31 'qssp.2010': 'fomosto_qssp2010', 

32 'qssp.2017': 'fomosto_qssp2017', 

33 'qssp.2020': 'fomosto_qssp2020', 

34 'qssp.ppeg2017': 'fomosto_qsspppeg2017', 

35} 

36 

37 

38def have_backend(version=None): 

39 avail = set() 

40 for version, exe in program_bins.items(): 

41 try: 

42 p = Popen([exe], stdout=PIPE, stderr=PIPE, stdin=PIPE) 

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

44 avail.add(version) 

45 

46 except OSError: 

47 pass 

48 

49 return avail 

50 

51 

52qssp_components = { 

53 1: 'ae an az gr sd te tn ue un uz ve vn vz'.split(), 

54 2: 'ar at ap gr sd tt tp ur ut up vr vt vp'.split(), 

55 3: '_disp_e _disp_n _disp_z'.split(), 

56 4: '_gravitation_e _gravitation_n _gravitation_z ' 

57 '_acce_e _acce_n _acce_z'.split(), 

58 5: '_rota_e _rota_n _rota_z'.split() 

59} 

60 

61 

62def str_float_vals(vals): 

63 return ' '.join(['%.6e' % val for val in vals]) 

64 

65 

66def cake_model_to_config(mod): 

67 k = 1000. 

68 srows = [] 

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

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

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

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

73 

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

75 

76 

77class QSSPSource(Object): 

78 lat = Float.T(default=0.0) 

79 lon = Float.T(default=0.0) 

80 depth = Float.T(default=10.0) 

81 torigin = Float.T(default=0.0) 

82 trise = Float.T(default=1.0) 

83 

84 def string_for_config(self): 

85 return '%(lat)15e %(lon)15e %(depth)15e %(torigin)15e %(trise)15e' \ 

86 % self.__dict__ 

87 

88 

89class QSSPSourceMT(QSSPSource): 

90 munit = Float.T(default=1.0) 

91 mrr = Float.T(default=1.0) 

92 mtt = Float.T(default=1.0) 

93 mpp = Float.T(default=1.0) 

94 mrt = Float.T(default=0.0) 

95 mrp = Float.T(default=0.0) 

96 mtp = Float.T(default=0.0) 

97 

98 def string_for_config(self): 

99 return '%(munit)15e %(mrr)15e %(mtt)15e %(mpp)15e ' \ 

100 '%(mrt)15e %(mrp)15e %(mtp)15e ' \ 

101 % self.__dict__ + QSSPSource.string_for_config(self) 

102 

103 

104class QSSPSourceDC(QSSPSource): 

105 moment = Float.T(default=1.0e9) 

106 strike = Float.T(default=0.0) 

107 dip = Float.T(default=90.0) 

108 rake = Float.T(default=0.0) 

109 

110 def string_for_config(self): 

111 return '%(moment)15e %(strike)15e %(dip)15e %(rake)15e ' \ 

112 % self.__dict__ + QSSPSource.string_for_config(self) 

113 

114 

115class QSSPReceiver(Object): 

116 lat = Float.T(default=10.0) 

117 lon = Float.T(default=0.0) 

118 name = String.T(default='') 

119 tstart = Float.T(default=0.0) 

120 distance = Float.T(default=0.0) 

121 

122 def string_for_config(self): 

123 return "%(lat)15e %(lon)15e '%(name)s' %(tstart)e" % self.__dict__ 

124 

125 

126class QSSPGreen(Object): 

127 depth = Float.T(default=10.0) 

128 filename = String.T(default='GF_10km') 

129 calculate = Bool.T(default=True) 

130 

131 def string_for_config(self): 

132 return "%(depth)15e '%(filename)s' %(calculate)i" % self.__dict__ 

133 

134 

135class QSSPConfig(Object): 

136 qssp_version = String.T(default='2010beta') 

137 time_region = Tuple.T(2, gf.Timing.T(), default=( 

138 gf.Timing('-10'), gf.Timing('+890'))) 

139 

140 frequency_max = Float.T(optional=True) 

141 slowness_max = Float.T(default=0.4) 

142 antialiasing_factor = Float.T(default=0.1) 

143 

144 # only available in 2017: 

145 switch_turning_point_filter = Int.T(default=0) 

146 max_pene_d1 = Float.T(default=2891.5) 

147 max_pene_d2 = Float.T(default=6371.0) 

148 earth_radius = Float.T(default=6371.0) 

149 switch_free_surf_reflection = Int.T(default=1) 

150 

151 lowpass_order = Int.T(default=0, optional=True) 

152 lowpass_corner = Float.T(default=1.0, optional=True) 

153 

154 bandpass_order = Int.T(default=0, optional=True) 

155 bandpass_corner_low = Float.T(default=1.0, optional=True) 

156 bandpass_corner_high = Float.T(default=1.0, optional=True) 

157 

158 output_slowness_min = Float.T(default=0.0, optional=True) 

159 output_slowness_max = Float.T(optional=True) 

160 

161 spheroidal_modes = Bool.T(default=True) 

162 toroidal_modes = Bool.T(default=True) 

163 

164 # only available in 2010beta: 

165 cutoff_harmonic_degree_sd = Int.T(optional=True, default=0) 

166 

167 cutoff_harmonic_degree_min = Int.T(default=0) 

168 cutoff_harmonic_degree_max = Int.T(default=25000) 

169 

170 crit_frequency_sge = Float.T(default=0.0) 

171 crit_harmonic_degree_sge = Int.T(default=0) 

172 

173 include_physical_dispersion = Bool.T(default=False) 

174 

175 source_patch_radius = Float.T(default=0.0) 

176 

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

178 

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

180 relevel_with_fade_in = Bool.T(default=False) 

181 nonzero_fade_in = Bool.T(default=False) 

182 nonzero_fade_out = Bool.T(default=False) 

183 

184 def items(self): 

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

186 

187 

188class QSSPConfigFull(QSSPConfig): 

189 time_window = Float.T(default=900.0) 

190 

191 receiver_depth = Float.T(default=0.0) 

192 sampling_interval = Float.T(default=5.0) 

193 

194 stored_quantity = String.T(default='displacement') 

195 output_filename = String.T(default='receivers') 

196 output_format = Int.T(default=1) 

197 output_time_window = Float.T(optional=True) 

198 

199 gf_directory = String.T(default='qssp_green') 

200 greens_functions = List.T(QSSPGreen.T()) 

201 

202 sources = List.T(QSSPSource.T()) 

203 receivers = List.T(QSSPReceiver.T()) 

204 

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

206 

207 @staticmethod 

208 def example(): 

209 conf = QSSPConfigFull() 

210 conf.sources.append(QSSPSourceMT()) 

211 lats = [20.] 

212 conf.receivers.extend(QSSPReceiver(lat=lat) for lat in lats) 

213 conf.greens_functions.append(QSSPGreen()) 

214 return conf 

215 

216 @property 

217 def components(self): 

218 if self.qssp_version in ('2017', '2020'): 

219 if self.stored_quantity == 'rotation': 

220 fmt = 5 

221 else: 

222 fmt = 3 

223 elif self.qssp_version == 'ppeg2017': 

224 fmt = 4 

225 else: 

226 fmt = self.output_format 

227 

228 return qssp_components[fmt] 

229 

230 def get_output_filenames(self, rundir): 

231 if self.qssp_version in ('2017', '2020', 'ppeg2017'): 

232 return [ 

233 pjoin(rundir, self.output_filename + c + '.dat') 

234 for c in self.components] 

235 else: 

236 return [ 

237 pjoin(rundir, self.output_filename + '.' + c) 

238 for c in self.components] 

239 

240 def ensure_gf_directory(self): 

241 util.ensuredir(self.gf_directory) 

242 

243 def string_for_config(self): 

244 

245 def aggregate(xx): 

246 return len(xx), '\n'.join(x.string_for_config() for x in xx) 

247 

248 assert len(self.greens_functions) > 0 

249 assert len(self.sources) > 0 

250 assert len(self.receivers) > 0 

251 

252 d = self.__dict__.copy() 

253 

254 if self.output_time_window is None: 

255 d['output_time_window'] = self.time_window 

256 

257 if self.output_slowness_max is None: 

258 d['output_slowness_max'] = self.slowness_max 

259 

260 if self.frequency_max is None: 

261 d['frequency_max'] = 0.5/self.sampling_interval 

262 

263 d['gf_directory'] = os.path.abspath(self.gf_directory) + '/' 

264 

265 d['n_receiver_lines'], d['receiver_lines'] = aggregate(self.receivers) 

266 d['n_source_lines'], d['source_lines'] = aggregate(self.sources) 

267 d['n_gf_lines'], d['gf_lines'] = aggregate(self.greens_functions) 

268 model_str, nlines = cake_model_to_config(self.earthmodel_1d) 

269 d['n_model_lines'] = nlines 

270 d['model_lines'] = model_str 

271 if self.stored_quantity == 'rotation': 

272 d['output_rotation'] = 1 

273 d['output_displacement'] = 0 

274 else: 

275 d['output_displacement'] = 1 

276 d['output_rotation'] = 0 

277 

278 if len(self.sources) == 0 or isinstance(self.sources[0], QSSPSourceMT): 

279 d['point_source_type'] = 1 

280 else: 

281 d['point_source_type'] = 2 

282 

283 if self.qssp_version == '2010beta': 

284 d['scutoff_doc'] = ''' 

285# (SH waves), and cutoff harmonic degree for static deformation 

286'''.strip() 

287 

288 d['scutoff'] = '%i' % self.cutoff_harmonic_degree_sd 

289 

290 d['sfilter_doc'] = ''' 

291# 3. selection of order of Butterworth low-pass filter (if <= 0, then no 

292# filtering), corner frequency (smaller than the cut-off frequency defined 

293# above) 

294'''.strip() 

295 

296 if self.bandpass_order != 0: 

297 raise QSSPError( 

298 'this version of qssp does not support bandpass ' 

299 'settings, use lowpass instead') 

300 

301 d['sfilter'] = '%i %f' % ( 

302 self.lowpass_order, 

303 self.lowpass_corner) 

304 

305 elif self.qssp_version in ('2010', '2017', '2020', 'ppeg2017'): 

306 d['scutoff_doc'] = ''' 

307# (SH waves), minimum and maximum cutoff harmonic degrees 

308# Note: if the near-field static displacement is desired, the minimum 

309# cutoff harmonic degree should not be smaller than, e.g., 2000. 

310'''.strip() 

311 

312 d['scutoff'] = '%i %i' % ( 

313 self.cutoff_harmonic_degree_min, 

314 self.cutoff_harmonic_degree_max) 

315 

316 d['sfilter_doc'] = ''' 

317# 3. selection of order of Butterworth bandpass filter (if <= 0, then no 

318# filtering), lower and upper corner frequencies (smaller than the cut-off 

319# frequency defined above) 

320'''.strip() 

321 

322 if self.lowpass_order != 0: 

323 raise QSSPError( 

324 'this version of qssp does not support lowpass settings, ' 

325 'use bandpass instead') 

326 

327 d['sfilter'] = '%i %f %f' % ( 

328 self.bandpass_order, 

329 self.bandpass_corner_low, 

330 self.bandpass_corner_high) 

331 if self.qssp_version in ('2017', '2020', 'ppeg2017'): 

332 template = '''# autogenerated QSSP input by qssp.py 

333# 

334# This is the input file of FORTRAN77 program "qssp2017" for calculating 

335# synthetic seismograms of a self-gravitating, spherically symmetric, 

336# isotropic and viscoelastic earth. 

337# 

338# by 

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

340# Helmholtz-Centre Potsdam 

341# GFZ German Reseach Centre for Geosciences 

342# Telegrafenberg, D-14473 Potsdam, Germany 

343# 

344# Last modified: Potsdam, October 2017 

345# 

346# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

347# If not specified, SI Unit System is used overall! 

348# 

349# Coordinate systems: 

350# spherical (r,t,p) with r = radial, 

351# t = co-latitude, 

352# p = east longitude. 

353# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

354# 

355# UNIFORM RECEIVER DEPTH 

356# ====================== 

357# 1. uniform receiver depth [km] 

358#------------------------------------------------------------------------------------------- 

359 %(receiver_depth)e 

360#------------------------------------------------------------------------------------------- 

361# 

362# SPACE-TIME SAMPLING PARAMETERS 

363# ========================= 

364# 1. time window [sec], sampling interval [sec] 

365# 2. max. frequency [Hz] of Green's functions 

366# 3. max. slowness [s/km] of Green's functions 

367# Note: if the near-field static displacement is desired, the maximum slowness should not 

368# be smaller than the S wave slowness in the receiver layer 

369# 4. anti-aliasing factor (> 0 & < 1), if it is <= 0 or >= 1/e (~ 0.4), then 

370# default value of 1/e is used (e.g., 0.1 = alias phases will be suppressed 

371# to 10%% of their original amplitude) 

372# 5. switch (1/0 = yes/no) of turning-point filter, the range (d1, d2) of max. penetration 

373# depth [km] (d1 is meaningless if it is smaller than the receiver/source depth, and 

374# d2 is meaningless if it is equal to or larger than the earth radius) 

375# 

376# Note: The turning-point filter (Line 5) works only for the extended QSSP code (e.g., 

377# qssp2016). if this filter is selected, all phases with the turning point 

378# shallower than d1 or deeper than d2 will be filtered. 

379# 

380# 6. Earth radius [km], switch of free-surface-reflection filter (1/0 = with/without free 

381# surface reflection) 

382# 

383# Note: The free-surface-reflection filter (Line 6) works only for the extended QSSP 

384# code (e.g., qssp2016). if this filter is selected, all phases with the turning 

385# point shallower than d1 or deeper than d2 will be filtered. 

386#------------------------------------------------------------------------------------------- 

387 %(time_window)e %(sampling_interval)e 

388 %(frequency_max)e 

389 %(slowness_max)e 

390 %(antialiasing_factor)e 

391 %(switch_turning_point_filter)i %(max_pene_d1)e %(max_pene_d2)e 

392 %(earth_radius)e %(switch_free_surf_reflection)i 

393#------------------------------------------------------------------------------------------- 

394# 

395# SELF-GRAVITATING EFFECT 

396# ======================= 

397# 1. the critical frequency [Hz] and the critical harmonic degree, below which 

398# the self-gravitating effect should be included 

399#------------------------------------------------------------------------------------------- 

400 %(crit_frequency_sge)e %(crit_harmonic_degree_sge)i 

401#------------------------------------------------------------------------------------------- 

402# 

403# WAVE TYPES 

404# ========== 

405# 1. selection (1/0 = yes/no) of speroidal modes (P-SV waves), selection of toroidal modes 

406%(scutoff_doc)s 

407#------------------------------------------------------------------------------------------- 

408 %(spheroidal_modes)i %(toroidal_modes)i %(scutoff)s 

409#------------------------------------------------------------------------------------------- 

410# GREEN'S FUNCTION FILES 

411# ====================== 

412# 1. number of discrete source depths, estimated radius of each source patch [km] and 

413# directory for Green's functions 

414# 2. list of the source depths [km], the respective file names of the Green's 

415# functions (spectra) and the switch number (0/1) (0 = do not calculate 

416# this Green's function because it exists already, 1 = calculate or update 

417# this Green's function. Note: update is required if any of the above 

418# parameters is changed) 

419#------------------------------------------------------------------------------------------- 

420 %(n_gf_lines)i %(source_patch_radius)e '%(gf_directory)s' 

421 %(gf_lines)s 

422#-------------------------------------------------------------------------------------------------------- 

423# 

424# MULTI-EVENT SOURCE PARAMETERS 

425# ============================= 

426# 1. number of discrete point sources and selection of the source data format 

427# (1, 2 or 3) 

428# 2. list of the multi-event sources 

429# 

430# Format 1 (full moment tensor): 

431# Unit Mrr Mtt Mpp Mrt Mrp Mtp Lat Lon Depth T_origin T_rise 

432# [Nm] [deg] [deg] [km] [sec] [sec] 

433# 

434# Format 2 (double couple): 

435# Unit Strike Dip Rake Lat Lon Depth T_origin T_rise 

436# [Nm] [deg] [deg] [deg] [deg] [deg] [km] [sec] [sec] 

437# 

438# Format 3 (single force): 

439# Unit Feast Fnorth Fvertical Lat Lon Depth T_origin T_rise 

440# [N] [deg] [deg] [km] [sec] [sec] 

441# 

442# Note: for each point source, the default moment (force) rate time function is used, defined by a 

443# squared half-period (T_rise) sinusoid starting at T_origin. 

444#----------------------------------------------------------------------------------- 

445 %(n_source_lines)i %(point_source_type)i 

446%(source_lines)s 

447#-------------------------------------------------------------------------------------------------------- 

448# 

449# RECEIVER PARAMETERS 

450# =================== 

451# 1. select output observables (1/0 = yes/no) 

452# Note: the gravity change defined here is space based, i.e., the effect due to free-air 

453# gradient and inertial are not included. the vertical component is positve upwards. 

454# 2. output file name 

455# 3. output time window [sec] (<= Green's function time window) 

456%(sfilter_doc)s 

457# 5. lower and upper slowness cut-off [s/km] (slowness band-pass filter) 

458# 6. number of receiver 

459# 7. list of the station parameters 

460# Format: 

461# Lat Lon Name Time_reduction 

462# [deg] [deg] [sec] 

463# (Note: Time_reduction = start time of the time window) 

464#--------------------------------------------------------------------------------------------------------------- 

465# disp | velo | acce | strain | strain_rate | stress | stress_rate | rotation | rot_rate | gravitation | gravity 

466#--------------------------------------------------------------------------------------------------------------- 

467# 1 1 1 1 1 1 1 1 1 1 1 

468 %(output_displacement)i 0 0 0 0 0 0 %(output_rotation)i 0 0 0 

469 

470 '%(output_filename)s' 

471 %(output_time_window)e 

472 %(sfilter)s 

473 %(output_slowness_min)e %(output_slowness_max)e 

474 %(n_receiver_lines)i 

475%(receiver_lines)s 

476#------------------------------------------------------------------------------------------- 

477# 

478# LAYERED EARTH MODEL (IASP91) 

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

480# 1. number of data lines of the layered model and selection for including 

481# the physical dispersion according Kamamori & Anderson (1977) 

482#------------------------------------------------------------------------------------------- 

483 %(n_model_lines)i %(include_physical_dispersion)i 

484#-------------------------------------------------------------------------------------------------------- 

485# 

486# MODEL PARAMETERS 

487# ================ 

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

489#------------------------------------------------------------------------------------------- 

490%(model_lines)s 

491#---------------------------------end of all inputs----------------------------------------- 

492''' # noqa 

493 else: 

494 

495 template = '''# autogenerated QSSP input by qssp.py 

496# 

497# This is the input file of FORTRAN77 program "qssp2010" for calculating 

498# synthetic seismograms of a self-gravitating, spherically symmetric, 

499# isotropic and viscoelastic earth. 

500# 

501# by 

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

503# Helmholtz-Centre Potsdam 

504# GFZ German Reseach Centre for Geosciences 

505# Telegrafenberg, D-14473 Potsdam, Germany 

506# 

507# Last modified: Potsdam, July, 2010 

508# 

509# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

510# If not specified, SI Unit System is used overall! 

511# 

512# Coordinate systems: 

513# spherical (r,t,p) with r = radial, 

514# t = co-latitude, 

515# p = east longitude. 

516# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

517# 

518# UNIFORM RECEIVER DEPTH 

519# ====================== 

520# 1. uniform receiver depth [km] 

521#------------------------------------------------------------------------------------------- 

522 %(receiver_depth)e 

523#------------------------------------------------------------------------------------------- 

524# 

525# TIME (FREQUENCY) SAMPLING 

526# ========================= 

527# 1. time window [sec], sampling interval [sec] 

528# 2. max. frequency [Hz] of Green's functions 

529# 3. max. slowness [s/km] of Green's functions 

530# Note: if the near-field static displacement is desired, the maximum slowness should not 

531# be smaller than the S wave slowness in the receiver layer 

532# 4. anti-aliasing factor (> 0 & < 1), if it is <= 0 or >= 1/e (~ 0.4), then 

533# default value of 1/e is used (e.g., 0.1 = alias phases will be suppressed 

534# to 10%% of their original amplitude) 

535# 

536# Note: The computation effort increases linearly the time window and 

537# quadratically with the cut-off frequency. 

538#------------------------------------------------------------------------------------------- 

539 %(time_window)e %(sampling_interval)e 

540 %(frequency_max)e 

541 %(slowness_max)e 

542 %(antialiasing_factor)e 

543#------------------------------------------------------------------------------------------- 

544# 

545# SELF-GRAVITATING EFFECT 

546# ======================= 

547# 1. the critical frequency [Hz] and the critical harmonic degree, below which 

548# the self-gravitating effect should be included 

549#------------------------------------------------------------------------------------------- 

550 %(crit_frequency_sge)e %(crit_harmonic_degree_sge)i 

551#------------------------------------------------------------------------------------------- 

552# 

553# WAVE TYPES 

554# ========== 

555# 1. selection (1/0 = yes/no) of speroidal modes (P-SV waves), selection of toroidal modes 

556%(scutoff_doc)s 

557#------------------------------------------------------------------------------------------- 

558 %(spheroidal_modes)i %(toroidal_modes)i %(scutoff)s 

559#------------------------------------------------------------------------------------------- 

560# GREEN'S FUNCTION FILES 

561# ====================== 

562# 1. number of discrete source depths, estimated radius of each source patch [km] and 

563# directory for Green's functions 

564# 2. list of the source depths [km], the respective file names of the Green's 

565# functions (spectra) and the switch number (0/1) (0 = do not calculate 

566# this Green's function because it exists already, 1 = calculate or update 

567# this Green's function. Note: update is required if any of the above 

568# parameters is changed) 

569#------------------------------------------------------------------------------------------- 

570 %(n_gf_lines)i %(source_patch_radius)e '%(gf_directory)s' 

571 %(gf_lines)s 

572#------------------------------------------------------------------------------------------- 

573# 

574# MULTI-EVENT SOURCE PARAMETERS 

575# ============================= 

576# 1. number of discrete point sources and selection of the source data format 

577# (1 or 2) 

578# 2. list of the multi-event sources 

579# Format 1: 

580# M-Unit Mrr Mtt Mpp Mrt Mrp Mtp Lat Lon Depth T_origin T_rise 

581# [Nm] [deg] [deg] [km] [sec] [sec] 

582# Format 2: 

583# Moment Strike Dip Rake Lat Lon Depth T_origin T_rise 

584# [Nm] [deg] [deg] [deg] [deg] [deg] [km] [sec] [sec] 

585#------------------------------------------------------------------------------------------- 

586 %(n_source_lines)i %(point_source_type)i 

587%(source_lines)s 

588#------------------------------------------------------------------------------------------- 

589# 

590# RECEIVER PARAMETERS 

591# =================== 

592# 1. output file name and and selection of output format: 

593# 1 = cartesian: vertical(z)/north(n)/east(e); 

594# 2 = spherical: radial(r)/theta(t)/phi(p) 

595# (Note: if output format 2 is selected, the epicenter (T_origin = 0) 

596# 2. output time window [sec] (<= Green's function time window) 

597%(sfilter_doc)s 

598# 4. lower and upper slowness cut-off [s/km] (slowness band-pass filter) 

599# 5. number of receiver 

600# 6. list of the station parameters 

601# Format: 

602# Lat Lon Name Time_reduction 

603# [deg] [deg] [sec] 

604# (Note: Time_reduction = start time of the time window) 

605#------------------------------------------------------------------------------------------- 

606 '%(output_filename)s' %(output_format)i 

607 %(output_time_window)e 

608 %(sfilter)s 

609 %(output_slowness_min)e %(output_slowness_max)e 

610 %(n_receiver_lines)i 

611%(receiver_lines)s 

612#------------------------------------------------------------------------------------------- 

613# 

614# LAYERED EARTH MODEL (IASP91) 

615# ============================ 

616# 1. number of data lines of the layered model and selection for including 

617# the physical dispersion according Kamamori & Anderson (1977) 

618#------------------------------------------------------------------------------------------- 

619 %(n_model_lines)i %(include_physical_dispersion)i 

620#------------------------------------------------------------------------------------------- 

621# 

622# MULTILAYERED MODEL PARAMETERS (source site) 

623# =========================================== 

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

625#------------------------------------------------------------------------------------------- 

626%(model_lines)s 

627#---------------------------------end of all inputs----------------------------------------- 

628''' # noqa 

629 

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

631 

632 

633class QSSPError(gf.store.StoreError): 

634 pass 

635 

636 

637class Interrupted(gf.store.StoreError): 

638 def __str__(self): 

639 return 'Interrupted.' 

640 

641 

642class QSSPRunner(object): 

643 

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

645 

646 self.tempdir = mkdtemp(prefix='qssprun-', dir=tmp) 

647 self.keep_tmp = keep_tmp 

648 self.config = None 

649 

650 def run(self, config): 

651 self.config = config 

652 

653 input_fn = pjoin(self.tempdir, 'input') 

654 

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

656 input_str = config.string_for_config() 

657 logger.debug('===== begin qssp input =====\n' 

658 '%s===== end qssp input =====' % input_str.decode()) 

659 f.write(input_str) 

660 

661 program = program_bins['qssp.%s' % config.qssp_version] 

662 

663 old_wd = os.getcwd() 

664 

665 os.chdir(self.tempdir) 

666 

667 interrupted = [] 

668 

669 def signal_handler(signum, frame): 

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

671 interrupted.append(True) 

672 

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

674 try: 

675 try: 

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

677 except OSError: 

678 os.chdir(old_wd) 

679 raise QSSPError( 

680 '''could not start qssp executable: "%s" 

681Available fomosto backends and download links to the modelling codes are listed 

682on 

683 

684 https://pyrocko.org/docs/current/apps/fomosto/backends.html 

685 

686''' % program) 

687 

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

689 

690 finally: 

691 signal.signal(signal.SIGINT, original) 

692 

693 if interrupted: 

694 raise KeyboardInterrupt() 

695 

696 logger.debug('===== begin qssp output =====\n' 

697 '%s===== end qssp output =====' % output_str.decode()) 

698 

699 errmess = [] 

700 if proc.returncode != 0: 

701 errmess.append( 

702 'qssp had a non-zero exit state: %i' % proc.returncode) 

703 if error_str: 

704 

705 logger.warning( 

706 'qssp emitted something via stderr: \n\n%s' 

707 % error_str.decode()) 

708 

709 # errmess.append('qssp emitted something via stderr') 

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

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

712 

713 if errmess: 

714 os.chdir(old_wd) 

715 raise QSSPError(''' 

716===== begin qssp input ===== 

717%s===== end qssp input ===== 

718===== begin qssp output ===== 

719%s===== end qssp output ===== 

720===== begin qssp error ===== 

721%s===== end qssp error ===== 

722%s 

723qssp has been invoked as "%s"'''.lstrip() % ( 

724 input_str.decode(), 

725 output_str.decode(), 

726 error_str.decode(), 

727 '\n'.join(errmess), 

728 program)) 

729 

730 self.qssp_output = output_str 

731 self.qssp_error = error_str 

732 

733 os.chdir(old_wd) 

734 

735 def get_traces(self): 

736 

737 fns = self.config.get_output_filenames(self.tempdir) 

738 traces = {} 

739 for comp, fn in zip(self.config.components, fns): 

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

741 nsamples, ntraces = data.shape 

742 ntraces -= 1 

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

744 toffset = data[0, 0] 

745 for itrace in range(ntraces): 

746 rec = self.config.receivers[itrace] 

747 tmin = rec.tstart + toffset 

748 tr = trace.Trace( 

749 '', '%04i' % itrace, '', comp, 

750 tmin=tmin, deltat=deltat, ydata=data[:, itrace+1], 

751 meta=dict(distance=rec.distance)) 

752 

753 traces[itrace, comp] = tr 

754 

755 if self.config.qssp_version == 'ppeg2017': 

756 for itrace in range(ntraces): 

757 for c in 'nez': 

758 tr_accel = traces[itrace, '_acce_' + c] 

759 tr_gravi = traces[itrace, '_gravitation_' + c] 

760 tr_ag = tr_accel.copy() 

761 tr_ag.ydata -= tr_gravi.ydata 

762 tr_ag.set_codes(channel='ag_' + c) 

763 tr_ag.meta = tr_accel.meta 

764 traces[itrace, 'ag_' + c] = tr_ag 

765 

766 traces = list(traces.values()) 

767 traces.sort(key=lambda tr: tr.nslc_id) 

768 return traces 

769 

770 def __del__(self): 

771 if self.tempdir: 

772 if not self.keep_tmp: 

773 shutil.rmtree(self.tempdir) 

774 self.tempdir = None 

775 else: 

776 logger.warning( 

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

778 

779 

780class QSSPGFBuilder(gf.builder.Builder): 

781 nsteps = 2 

782 

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

784 force=False): 

785 

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

787 baseconf = self.store.get_extra('qssp') 

788 if baseconf.qssp_version in ('2017', '2020'): 

789 if self.store.config.stored_quantity == 'rotation': 

790 self.gfmapping = [ 

791 (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)), 

792 {'_rota_n': (0, -1), '_rota_e': (3, -1), 

793 '_rota_z': (5, -1)}), 

794 (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)), 

795 {'_rota_n': (1, -1), '_rota_e': (4, -1), 

796 '_rota_z': (6, -1)}), 

797 (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)), 

798 {'_rota_n': (2, -1), '_rota_z': (7, -1)}), 

799 (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)), 

800 {'_rota_n': (8, -1), '_rota_z': (9, -1)}), 

801 ] 

802 else: 

803 self.gfmapping = [ 

804 (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)), 

805 {'_disp_n': (0, -1), '_disp_e': (3, -1), 

806 '_disp_z': (5, -1)}), 

807 (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)), 

808 {'_disp_n': (1, -1), '_disp_e': (4, -1), 

809 '_disp_z': (6, -1)}), 

810 (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)), 

811 {'_disp_n': (2, -1), '_disp_z': (7, -1)}), 

812 (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)), 

813 {'_disp_n': (8, -1), '_disp_z': (9, -1)}), 

814 ] 

815 

816 elif baseconf.qssp_version == 'ppeg2017': 

817 self.gfmapping = [ 

818 (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)), 

819 {'ag_n': (0, -1), 'ag_e': (3, -1), 'ag_z': (5, -1)}), 

820 (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)), 

821 {'ag_n': (1, -1), 'ag_e': (4, -1), 'ag_z': (6, -1)}), 

822 (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)), 

823 {'ag_n': (2, -1), 'ag_z': (7, -1)}), 

824 (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)), 

825 {'ag_n': (8, -1), 'ag_z': (9, -1)}), 

826 ] 

827 else: 

828 self.gfmapping = [ 

829 (MomentTensor(m=symmat6(1, 0, 0, 1, 0, 0)), 

830 {'un': (0, -1), 'ue': (3, -1), 'uz': (5, -1)}), 

831 (MomentTensor(m=symmat6(0, 0, 0, 0, 1, 1)), 

832 {'un': (1, -1), 'ue': (4, -1), 'uz': (6, -1)}), 

833 (MomentTensor(m=symmat6(0, 0, 1, 0, 0, 0)), 

834 {'un': (2, -1), 'uz': (7, -1)}), 

835 (MomentTensor(m=symmat6(0, 1, 0, 0, 0, 0)), 

836 {'un': (8, -1), 'uz': (9, -1)}), 

837 ] 

838 

839 if step == 0: 

840 block_size = (1, 1, self.store.config.ndistances)