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) 

841 else: 

842 if block_size is None: 

843 block_size = (1, 1, 51) 

844 

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

846 block_size = block_size[1:] 

847 

848 gf.builder.Builder.__init__( 

849 self, self.store.config, step, block_size=block_size, force=force) 

850 

851 conf = QSSPConfigFull(**baseconf.items()) 

852 conf.gf_directory = pjoin(store_dir, 'qssp_green') 

853 conf.earthmodel_1d = self.store.config.earthmodel_1d 

854 deltat = self.store.config.deltat 

855 if 'time_window' not in shared: 

856 d = self.store.make_timing_params( 

857 conf.time_region[0], conf.time_region[1], 

858 force=force) 

859 

860 tmax = math.ceil(d['tmax'] / deltat) * deltat 

861 tmin = math.floor(d['tmin'] / deltat) * deltat 

862 

863 shared['time_window'] = tmax - tmin 

864 shared['tstart'] = tmin 

865 

866 self.tstart = shared['tstart'] 

867 conf.time_window = shared['time_window'] 

868 

869 self.tmp = tmp 

870 if self.tmp is not None: 

871 util.ensuredir(self.tmp) 

872 

873 util.ensuredir(conf.gf_directory) 

874 

875 self.qssp_config = conf 

876 

877 def cleanup(self): 

878 self.store.close() 

879 

880 def work_block(self, iblock): 

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

882 (sz, firstx), (sz, lastx), (ns, nx) = \ 

883 self.get_block_extents(iblock) 

884 

885 rz = self.store.config.receiver_depth 

886 else: 

887 (rz, sz, firstx), (rz, sz, lastx), (nr, ns, nx) = \ 

888 self.get_block_extents(iblock) 

889 

890 gf_filename = 'GF_%gkm_%gkm' % (sz/km, rz/km) 

891 

892 conf = copy.deepcopy(self.qssp_config) 

893 

894 gf_path = os.path.join(conf.gf_directory, '?_' + gf_filename) 

895 

896 if self.step == 0 and len(glob.glob(gf_path)) > 0: 

897 logger.info( 

898 'Skipping step %i / %i, block %i / %i (GF already exists)' 

899 % (self.step+1, self.nsteps, iblock+1, self.nblocks)) 

900 

901 return 

902 

903 logger.info( 

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

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

906 

907 tbeg = time.time() 

908 

909 runner = QSSPRunner(tmp=self.tmp) 

910 

911 conf.receiver_depth = rz/km 

912 conf.stored_quantity = self.store.config.stored_quantity 

913 conf.sampling_interval = 1.0 / self.gf_config.sample_rate 

914 dx = self.gf_config.distance_delta 

915 

916 if self.step == 0: 

917 distances = [firstx] 

918 else: 

919 distances = num.linspace(firstx, firstx + (nx-1)*dx, nx) 

920 

921 conf.receivers = [ 

922 QSSPReceiver( 

923 lat=90-d*cake.m2d, 

924 lon=180., 

925 tstart=self.tstart, 

926 distance=d) 

927 

928 for d in distances] 

929 

930 if self.step == 0: 

931 gf_filename = 'TEMP' + gf_filename[2:] 

932 

933 gfs = [QSSPGreen( 

934 filename=gf_filename, 

935 depth=sz/km, 

936 calculate=(self.step == 0))] 

937 

938 conf.greens_functions = gfs 

939 

940 trise = 0.001*conf.sampling_interval # make it short (delta impulse) 

941 

942 if self.step == 0: 

943 conf.sources = [QSSPSourceMT( 

944 lat=90-0.001*dx*cake.m2d, 

945 lon=0.0, 

946 trise=trise, 

947 torigin=0.0)] 

948 

949 runner.run(conf) 

950 gf_path = os.path.join(conf.gf_directory, '?_' + gf_filename) 

951 for s in glob.glob(gf_path): 

952 d = s.replace('TEMP_', 'GF_') 

953 os.rename(s, d) 

954 

955 else: 

956 for mt, gfmap in self.gfmapping[ 

957 :[3, 4][self.gf_config.ncomponents == 10]]: 

958 m = mt.m_up_south_east() 

959 

960 conf.sources = [QSSPSourceMT( 

961 lat=90-0.001*dx*cake.m2d, 

962 lon=0.0, 

963 mrr=m[0, 0], mtt=m[1, 1], mpp=m[2, 2], 

964 mrt=m[0, 1], mrp=m[0, 2], mtp=m[1, 2], 

965 trise=trise, 

966 torigin=0.0)] 

967 

968 runner.run(conf) 

969 

970 rawtraces = runner.get_traces() 

971 

972 interrupted = [] 

973 

974 def signal_handler(signum, frame): 

975 interrupted.append(True) 

976 

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

978 self.store.lock() 

979 duplicate_inserts = 0 

980 try: 

981 for itr, tr in enumerate(rawtraces): 

982 if tr.channel in gfmap: 

983 

984 x = tr.meta['distance'] 

985 ig, factor = gfmap[tr.channel] 

986 

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

988 args = (sz, x, ig) 

989 else: 

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

991 

992 if conf.cut: 

993 tmin = self.store.t(conf.cut[0], args[:-1]) 

994 tmax = self.store.t(conf.cut[1], args[:-1]) 

995 if None in (tmin, tmax): 

996 continue 

997 

998 tr.chop(tmin, tmax) 

999 

1000 tmin = tr.tmin 

1001 tmax = tr.tmax 

1002 

1003 if conf.fade: 

1004 ta, tb, tc, td = [ 

1005 self.store.t(v, args[:-1]) 

1006 for v in conf.fade] 

1007 

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

1009 continue 

1010 

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

1012 raise QSSPError( 

1013 'invalid fade configuration') 

1014 

1015 t = tr.get_xdata() 

1016 fin = num.interp(t, [ta, tb], [0., 1.]) 

1017 fout = num.interp(t, [tc, td], [1., 0.]) 

1018 anti_fin = 1. - fin 

1019 anti_fout = 1. - fout 

1020 

1021 y = tr.ydata 

1022 

1023 sum_anti_fin = num.sum(anti_fin) 

1024 sum_anti_fout = num.sum(anti_fout) 

1025 

1026 if conf.nonzero_fade_in \ 

1027 and sum_anti_fin != 0.0: 

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

1029 else: 

1030 yin = 0.0 

1031 

1032 if conf.nonzero_fade_out \ 

1033 and sum_anti_fout != 0.0: 

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

1035 else: 

1036 yout = 0.0 

1037 

1038 y2 = anti_fin*yin + fin*fout*y + anti_fout*yout 

1039 

1040 if conf.relevel_with_fade_in: 

1041 y2 -= yin 

1042 

1043 tr.set_ydata(y2) 

1044 

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

1046 gf_tr.data *= factor 

1047 

1048 try: 

1049 self.store.put(args, gf_tr) 

1050 except gf.store.DuplicateInsert: 

1051 duplicate_inserts += 1 

1052 

1053 finally: 

1054 if duplicate_inserts: 

1055 logger.warning( 

1056 '%i insertions skipped (duplicates)' 

1057 % duplicate_inserts) 

1058 

1059 self.store.unlock() 

1060 signal.signal(signal.SIGINT, original) 

1061 

1062 if interrupted: 

1063 raise KeyboardInterrupt() 

1064 

1065 tend = time.time() 

1066 logger.info( 

1067 'Done with step %i / %i, block %i / %i, wallclock time: %.0f s' % ( 

1068 self.step+1, self.nsteps, iblock+1, self.nblocks, tend-tbeg)) 

1069 

1070 

1071km = 1000. 

1072 

1073 

1074def init(store_dir, variant, config_params=None): 

1075 if variant is None: 

1076 variant = '2010beta' 

1077 

1078 if ('qssp.' + variant) not in program_bins: 

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

1080 

1081 qssp = QSSPConfig(qssp_version=variant) 

1082 if variant != 'ppeg2017': 

1083 qssp.time_region = ( 

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

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

1086 

1087 qssp.cut = ( 

1088 gf.Timing('begin-50'), 

1089 gf.Timing('end+100')) 

1090 

1091 else: # variant == 'ppeg2017': 

1092 qssp.frequency_max = 0.5 

1093 qssp.time_region = [ 

1094 gf.Timing('-100'), gf.Timing('{stored:begin}+100')] 

1095 qssp.cut = [ 

1096 gf.Timing('-100'), gf.Timing('{stored:begin}+100')] 

1097 qssp.antialiasing_factor = 1.0e-10 

1098 qssp.toroidal_modes = False 

1099 qssp.cutoff_harmonic_degree_min = 2500 

1100 qssp.cutoff_harmonic_degree_max = 2500 

1101 qssp.crit_frequency_sge = 5.0 

1102 qssp.crit_harmonic_degree_sge = 50000 

1103 qssp.source_patch_radius = 10.0 

1104 qssp.bandpass_order = 6 

1105 qssp.bandpass_corner_low = 0.0 

1106 qssp.bandpass_corner_high = 0.125 

1107 

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

1109 if variant == 'ppeg2017': 

1110 quantity = 'acceleration' 

1111 else: 

1112 quantity = None 

1113 

1114 if variant == 'ppeg2017': 

1115 sample_rate = 4.0 

1116 else: 

1117 sample_rate = 0.2 

1118 

1119 d = dict( 

1120 id=store_id, 

1121 ncomponents=10, 

1122 component_scheme='elastic10', 

1123 stored_quantity=quantity, 

1124 sample_rate=sample_rate, 

1125 receiver_depth=0*km, 

1126 source_depth_min=10*km, 

1127 source_depth_max=20*km, 

1128 source_depth_delta=10*km, 

1129 distance_min=100*km, 

1130 distance_max=1000*km, 

1131 distance_delta=10*km, 

1132 earthmodel_1d=cake.load_model(), 

1133 modelling_code_id='qssp', 

1134 tabulated_phases=[ 

1135 gf.meta.TPDef( 

1136 id='begin', 

1137 definition='p,P,p\\,P\\,Pv_(cmb)p'), 

1138 gf.meta.TPDef( 

1139 id='end', 

1140 definition='2.5'), 

1141 gf.meta.TPDef( 

1142 id='P', 

1143 definition='!P'), 

1144 gf.meta.TPDef( 

1145 id='S', 

1146 definition='!S'), 

1147 gf.meta.TPDef( 

1148 id='p', 

1149 definition='!p'), 

1150 gf.meta.TPDef( 

1151 id='s', 

1152 definition='!s')]) 

1153 

1154 if config_params is not None: 

1155 d.update(config_params) 

1156 

1157 config = gf.meta.ConfigTypeA(**d) 

1158 

1159 config.validate() 

1160 return gf.store.Store.create_editables( 

1161 store_dir, 

1162 config=config, 

1163 extra={'qssp': qssp}) 

1164 

1165 

1166def build( 

1167 store_dir, 

1168 force=False, 

1169 nworkers=None, 

1170 continue_=False, 

1171 step=None, 

1172 iblock=None): 

1173 

1174 return QSSPGFBuilder.build( 

1175 store_dir, force=force, nworkers=nworkers, continue_=continue_, 

1176 step=step, iblock=iblock) 

1177 

1178 

1179def get_conf( 

1180 store_dir, 

1181 force=False, 

1182 nworkers=None, 

1183 continue_=False, 

1184 step=None, 

1185 iblock=None): 

1186 

1187 return QSSPGFBuilder.get_conf()