1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5from __future__ import absolute_import, division 

6 

7import numpy as num 

8import logging 

9import os 

10import shutil 

11import glob 

12import copy 

13import signal 

14import math 

15import time 

16 

17from tempfile import mkdtemp 

18from subprocess import Popen, PIPE 

19from os.path import join as pjoin 

20 

21from pyrocko import trace, util, cake, gf 

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

23from pyrocko.moment_tensor import MomentTensor, symmat6 

24 

25guts_prefix = 'pf' 

26 

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

28 

29# how to call the programs 

30program_bins = { 

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

32 'qssp.2010': 'fomosto_qssp2010', 

33 'qssp.2017': 'fomosto_qssp2017', 

34 'qssp.2020': 'fomosto_qssp2020', 

35 'qssp.ppeg2017': 'fomosto_qsspppeg2017', 

36} 

37 

38 

39def have_backend(version=None): 

40 avail = set() 

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

42 try: 

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

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

45 avail.add(version) 

46 

47 except OSError: 

48 pass 

49 

50 return avail 

51 

52 

53qssp_components = { 

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

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

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

57 4: '_gravitation_e _gravitation_n _gravitation_z ' 

58 '_acce_e _acce_n _acce_z'.split(), 

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

60} 

61 

62 

63def str_float_vals(vals): 

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

65 

66 

67def cake_model_to_config(mod): 

68 k = 1000. 

69 srows = [] 

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

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

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

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

74 

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

76 

77 

78class QSSPSource(Object): 

79 lat = Float.T(default=0.0) 

80 lon = Float.T(default=0.0) 

81 depth = Float.T(default=10.0) 

82 torigin = Float.T(default=0.0) 

83 trise = Float.T(default=1.0) 

84 

85 def string_for_config(self): 

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

87 % self.__dict__ 

88 

89 

90class QSSPSourceMT(QSSPSource): 

91 munit = Float.T(default=1.0) 

92 mrr = Float.T(default=1.0) 

93 mtt = Float.T(default=1.0) 

94 mpp = Float.T(default=1.0) 

95 mrt = Float.T(default=0.0) 

96 mrp = Float.T(default=0.0) 

97 mtp = Float.T(default=0.0) 

98 

99 def string_for_config(self): 

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

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

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

103 

104 

105class QSSPSourceDC(QSSPSource): 

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

107 strike = Float.T(default=0.0) 

108 dip = Float.T(default=90.0) 

109 rake = Float.T(default=0.0) 

110 

111 def string_for_config(self): 

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

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

114 

115 

116class QSSPReceiver(Object): 

117 lat = Float.T(default=10.0) 

118 lon = Float.T(default=0.0) 

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

120 tstart = Float.T(default=0.0) 

121 distance = Float.T(default=0.0) 

122 

123 def string_for_config(self): 

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

125 

126 

127class QSSPGreen(Object): 

128 depth = Float.T(default=10.0) 

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

130 calculate = Bool.T(default=True) 

131 

132 def string_for_config(self): 

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

134 

135 

136class QSSPConfig(Object): 

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

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

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

140 

141 frequency_max = Float.T(optional=True) 

142 slowness_max = Float.T(default=0.4) 

143 antialiasing_factor = Float.T(default=0.1) 

144 

145 # only available in 2017: 

146 switch_turning_point_filter = Int.T(default=0) 

147 max_pene_d1 = Float.T(default=2891.5) 

148 max_pene_d2 = Float.T(default=6371.0) 

149 earth_radius = Float.T(default=6371.0) 

150 switch_free_surf_reflection = Int.T(default=1) 

151 

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

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

154 

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

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

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

158 

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

160 output_slowness_max = Float.T(optional=True) 

161 

162 spheroidal_modes = Bool.T(default=True) 

163 toroidal_modes = Bool.T(default=True) 

164 

165 # only available in 2010beta: 

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

167 

168 cutoff_harmonic_degree_min = Int.T(default=0) 

169 cutoff_harmonic_degree_max = Int.T(default=25000) 

170 

171 crit_frequency_sge = Float.T(default=0.0) 

172 crit_harmonic_degree_sge = Int.T(default=0) 

173 

174 include_physical_dispersion = Bool.T(default=False) 

175 

176 source_patch_radius = Float.T(default=0.0) 

177 

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

179 

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

181 relevel_with_fade_in = Bool.T(default=False) 

182 nonzero_fade_in = Bool.T(default=False) 

183 nonzero_fade_out = Bool.T(default=False) 

184 

185 def items(self): 

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

187 

188 

189class QSSPConfigFull(QSSPConfig): 

190 time_window = Float.T(default=900.0) 

191 

192 receiver_depth = Float.T(default=0.0) 

193 sampling_interval = Float.T(default=5.0) 

194 

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

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

197 output_format = Int.T(default=1) 

198 output_time_window = Float.T(optional=True) 

199 

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

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

202 

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

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

205 

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

207 

208 @staticmethod 

209 def example(): 

210 conf = QSSPConfigFull() 

211 conf.sources.append(QSSPSourceMT()) 

212 lats = [20.] 

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

214 conf.greens_functions.append(QSSPGreen()) 

215 return conf 

216 

217 @property 

218 def components(self): 

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

220 if self.stored_quantity == "rotation": 

221 fmt = 5 

222 else: 

223 fmt = 3 

224 elif self.qssp_version == 'ppeg2017': 

225 fmt = 4 

226 else: 

227 fmt = self.output_format 

228 

229 return qssp_components[fmt] 

230 

231 def get_output_filenames(self, rundir): 

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

233 return [ 

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

235 for c in self.components] 

236 else: 

237 return [ 

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

239 for c in self.components] 

240 

241 def ensure_gf_directory(self): 

242 util.ensuredir(self.gf_directory) 

243 

244 def string_for_config(self): 

245 

246 def aggregate(xx): 

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

248 

249 assert len(self.greens_functions) > 0 

250 assert len(self.sources) > 0 

251 assert len(self.receivers) > 0 

252 

253 d = self.__dict__.copy() 

254 

255 if self.output_time_window is None: 

256 d['output_time_window'] = self.time_window 

257 

258 if self.output_slowness_max is None: 

259 d['output_slowness_max'] = self.slowness_max 

260 

261 if self.frequency_max is None: 

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

263 

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

265 

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

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

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

269 model_str, nlines = cake_model_to_config(self.earthmodel_1d) 

270 d['n_model_lines'] = nlines 

271 d['model_lines'] = model_str 

272 if self.stored_quantity == "rotation": 

273 d['output_rotation'] = 1 

274 d['output_displacement'] = 0 

275 else: 

276 d['output_displacement'] = 1 

277 d['output_rotation'] = 0 

278 

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

280 d['point_source_type'] = 1 

281 else: 

282 d['point_source_type'] = 2 

283 

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

285 d['scutoff_doc'] = ''' 

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

287'''.strip() 

288 

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

290 

291 d['sfilter_doc'] = ''' 

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

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

294# above) 

295'''.strip() 

296 

297 if self.bandpass_order != 0: 

298 raise QSSPError( 

299 'this version of qssp does not support bandpass ' 

300 'settings, use lowpass instead') 

301 

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

303 self.lowpass_order, 

304 self.lowpass_corner) 

305 

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

307 d['scutoff_doc'] = ''' 

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

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

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

311'''.strip() 

312 

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

314 self.cutoff_harmonic_degree_min, 

315 self.cutoff_harmonic_degree_max) 

316 

317 d['sfilter_doc'] = ''' 

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

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

320# frequency defined above) 

321'''.strip() 

322 

323 if self.lowpass_order != 0: 

324 raise QSSPError( 

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

326 'use bandpass instead') 

327 

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

329 self.bandpass_order, 

330 self.bandpass_corner_low, 

331 self.bandpass_corner_high) 

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

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

334# 

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

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

337# isotropic and viscoelastic earth. 

338# 

339# by 

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

341# Helmholtz-Centre Potsdam 

342# GFZ German Reseach Centre for Geosciences 

343# Telegrafenberg, D-14473 Potsdam, Germany 

344# 

345# Last modified: Potsdam, October 2017 

346# 

347# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

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

349# 

350# Coordinate systems: 

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

352# t = co-latitude, 

353# p = east longitude. 

354# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

355# 

356# UNIFORM RECEIVER DEPTH 

357# ====================== 

358# 1. uniform receiver depth [km] 

359#------------------------------------------------------------------------------------------- 

360 %(receiver_depth)e 

361#------------------------------------------------------------------------------------------- 

362# 

363# SPACE-TIME SAMPLING PARAMETERS 

364# ========================= 

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

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

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

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

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

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

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

372# to 10%% of their original amplitude) 

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

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

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

376# 

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

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

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

380# 

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

382# surface reflection) 

383# 

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

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

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

387#------------------------------------------------------------------------------------------- 

388 %(time_window)e %(sampling_interval)e 

389 %(frequency_max)e 

390 %(slowness_max)e 

391 %(antialiasing_factor)e 

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

393 %(earth_radius)e %(switch_free_surf_reflection)i 

394#------------------------------------------------------------------------------------------- 

395# 

396# SELF-GRAVITATING EFFECT 

397# ======================= 

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

399# the self-gravitating effect should be included 

400#------------------------------------------------------------------------------------------- 

401 %(crit_frequency_sge)e %(crit_harmonic_degree_sge)i 

402#------------------------------------------------------------------------------------------- 

403# 

404# WAVE TYPES 

405# ========== 

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

407%(scutoff_doc)s 

408#------------------------------------------------------------------------------------------- 

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

410#------------------------------------------------------------------------------------------- 

411# GREEN'S FUNCTION FILES 

412# ====================== 

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

414# directory for Green's functions 

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

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

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

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

419# parameters is changed) 

420#------------------------------------------------------------------------------------------- 

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

422 %(gf_lines)s 

423#-------------------------------------------------------------------------------------------------------- 

424# 

425# MULTI-EVENT SOURCE PARAMETERS 

426# ============================= 

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

428# (1, 2 or 3) 

429# 2. list of the multi-event sources 

430# 

431# Format 1 (full moment tensor): 

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

433# [Nm] [deg] [deg] [km] [sec] [sec] 

434# 

435# Format 2 (double couple): 

436# Unit Strike Dip Rake Lat Lon Depth T_origin T_rise 

437# [Nm] [deg] [deg] [deg] [deg] [deg] [km] [sec] [sec] 

438# 

439# Format 3 (single force): 

440# Unit Feast Fnorth Fvertical Lat Lon Depth T_origin T_rise 

441# [N] [deg] [deg] [km] [sec] [sec] 

442# 

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

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

445#----------------------------------------------------------------------------------- 

446 %(n_source_lines)i %(point_source_type)i 

447%(source_lines)s 

448#-------------------------------------------------------------------------------------------------------- 

449# 

450# RECEIVER PARAMETERS 

451# =================== 

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

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

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

455# 2. output file name 

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

457%(sfilter_doc)s 

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

459# 6. number of receiver 

460# 7. list of the station parameters 

461# Format: 

462# Lat Lon Name Time_reduction 

463# [deg] [deg] [sec] 

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

465#--------------------------------------------------------------------------------------------------------------- 

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

467#--------------------------------------------------------------------------------------------------------------- 

468# 1 1 1 1 1 1 1 1 1 1 1 

469 %(output_displacement)i 0 0 0 0 0 0 %(output_rotation)i 0 0 0 

470 

471 '%(output_filename)s' 

472 %(output_time_window)e 

473 %(sfilter)s 

474 %(output_slowness_min)e %(output_slowness_max)e 

475 %(n_receiver_lines)i 

476%(receiver_lines)s 

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

478# 

479# LAYERED EARTH MODEL (IASP91) 

480# ============================ 

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

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

483#------------------------------------------------------------------------------------------- 

484 %(n_model_lines)i %(include_physical_dispersion)i 

485#-------------------------------------------------------------------------------------------------------- 

486# 

487# MODEL PARAMETERS 

488# ================ 

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

490#------------------------------------------------------------------------------------------- 

491%(model_lines)s 

492#---------------------------------end of all inputs----------------------------------------- 

493''' # noqa 

494 else: 

495 

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

497# 

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

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

500# isotropic and viscoelastic earth. 

501# 

502# by 

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

504# Helmholtz-Centre Potsdam 

505# GFZ German Reseach Centre for Geosciences 

506# Telegrafenberg, D-14473 Potsdam, Germany 

507# 

508# Last modified: Potsdam, July, 2010 

509# 

510# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

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

512# 

513# Coordinate systems: 

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

515# t = co-latitude, 

516# p = east longitude. 

517# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = 

518# 

519# UNIFORM RECEIVER DEPTH 

520# ====================== 

521# 1. uniform receiver depth [km] 

522#------------------------------------------------------------------------------------------- 

523 %(receiver_depth)e 

524#------------------------------------------------------------------------------------------- 

525# 

526# TIME (FREQUENCY) SAMPLING 

527# ========================= 

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

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

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

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

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

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

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

535# to 10%% of their original amplitude) 

536# 

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

538# quadratically with the cut-off frequency. 

539#------------------------------------------------------------------------------------------- 

540 %(time_window)e %(sampling_interval)e 

541 %(frequency_max)e 

542 %(slowness_max)e 

543 %(antialiasing_factor)e 

544#------------------------------------------------------------------------------------------- 

545# 

546# SELF-GRAVITATING EFFECT 

547# ======================= 

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

549# the self-gravitating effect should be included 

550#------------------------------------------------------------------------------------------- 

551 %(crit_frequency_sge)e %(crit_harmonic_degree_sge)i 

552#------------------------------------------------------------------------------------------- 

553# 

554# WAVE TYPES 

555# ========== 

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

557%(scutoff_doc)s 

558#------------------------------------------------------------------------------------------- 

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

560#------------------------------------------------------------------------------------------- 

561# GREEN'S FUNCTION FILES 

562# ====================== 

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

564# directory for Green's functions 

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

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

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

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

569# parameters is changed) 

570#------------------------------------------------------------------------------------------- 

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

572 %(gf_lines)s 

573#------------------------------------------------------------------------------------------- 

574# 

575# MULTI-EVENT SOURCE PARAMETERS 

576# ============================= 

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

578# (1 or 2) 

579# 2. list of the multi-event sources 

580# Format 1: 

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

582# [Nm] [deg] [deg] [km] [sec] [sec] 

583# Format 2: 

584# Moment Strike Dip Rake Lat Lon Depth T_origin T_rise 

585# [Nm] [deg] [deg] [deg] [deg] [deg] [km] [sec] [sec] 

586#------------------------------------------------------------------------------------------- 

587 %(n_source_lines)i %(point_source_type)i 

588%(source_lines)s 

589#------------------------------------------------------------------------------------------- 

590# 

591# RECEIVER PARAMETERS 

592# =================== 

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

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

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

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

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

598%(sfilter_doc)s 

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

600# 5. number of receiver 

601# 6. list of the station parameters 

602# Format: 

603# Lat Lon Name Time_reduction 

604# [deg] [deg] [sec] 

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

606#------------------------------------------------------------------------------------------- 

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

608 %(output_time_window)e 

609 %(sfilter)s 

610 %(output_slowness_min)e %(output_slowness_max)e 

611 %(n_receiver_lines)i 

612%(receiver_lines)s 

613#------------------------------------------------------------------------------------------- 

614# 

615# LAYERED EARTH MODEL (IASP91) 

616# ============================ 

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

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

619#------------------------------------------------------------------------------------------- 

620 %(n_model_lines)i %(include_physical_dispersion)i 

621#------------------------------------------------------------------------------------------- 

622# 

623# MULTILAYERED MODEL PARAMETERS (source site) 

624# =========================================== 

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

626#------------------------------------------------------------------------------------------- 

627%(model_lines)s 

628#---------------------------------end of all inputs----------------------------------------- 

629''' # noqa 

630 

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

632 

633 

634class QSSPError(gf.store.StoreError): 

635 pass 

636 

637 

638class Interrupted(gf.store.StoreError): 

639 def __str__(self): 

640 return 'Interrupted.' 

641 

642 

643class QSSPRunner(object): 

644 

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

646 

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

648 self.keep_tmp = keep_tmp 

649 self.config = None 

650 

651 def run(self, config): 

652 self.config = config 

653 

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

655 

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

657 input_str = config.string_for_config() 

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

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

660 f.write(input_str) 

661 

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

663 

664 old_wd = os.getcwd() 

665 

666 os.chdir(self.tempdir) 

667 

668 interrupted = [] 

669 

670 def signal_handler(signum, frame): 

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

672 interrupted.append(True) 

673 

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

675 try: 

676 try: 

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

678 except OSError: 

679 os.chdir(old_wd) 

680 raise QSSPError( 

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

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

683on 

684 

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

686 

687''' % program) 

688 

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

690 

691 finally: 

692 signal.signal(signal.SIGINT, original) 

693 

694 if interrupted: 

695 raise KeyboardInterrupt() 

696 

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

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

699 

700 errmess = [] 

701 if proc.returncode != 0: 

702 errmess.append( 

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

704 if error_str: 

705 

706 logger.warn( 

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

708 % error_str.decode()) 

709 

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

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

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

713 

714 if errmess: 

715 os.chdir(old_wd) 

716 raise QSSPError(''' 

717===== begin qssp input ===== 

718%s===== end qssp input ===== 

719===== begin qssp output ===== 

720%s===== end qssp output ===== 

721===== begin qssp error ===== 

722%s===== end qssp error ===== 

723%s 

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

725 input_str.decode(), 

726 output_str.decode(), 

727 error_str.decode(), 

728 '\n'.join(errmess), 

729 program)) 

730 

731 self.qssp_output = output_str 

732 self.qssp_error = error_str 

733 

734 os.chdir(old_wd) 

735 

736 def get_traces(self): 

737 

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

739 traces = {} 

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

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

742 nsamples, ntraces = data.shape 

743 ntraces -= 1 

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

745 toffset = data[0, 0] 

746 for itrace in range(ntraces): 

747 rec = self.config.receivers[itrace] 

748 tmin = rec.tstart + toffset 

749 tr = trace.Trace( 

750 '', '%04i' % itrace, '', comp, 

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

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

753 

754 traces[itrace, comp] = tr 

755 

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

757 for itrace in range(ntraces): 

758 for c in 'nez': 

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

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

761 tr_ag = tr_accel.copy() 

762 tr_ag.ydata -= tr_gravi.ydata 

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

764 tr_ag.meta = tr_accel.meta 

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

766 

767 traces = list(traces.values()) 

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

769 return traces 

770 

771 def __del__(self): 

772 if self.tempdir: 

773 if not self.keep_tmp: 

774 shutil.rmtree(self.tempdir) 

775 self.tempdir = None 

776 else: 

777 logger.warn( 

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

779 

780 

781class QSSPGFBuilder(gf.builder.Builder): 

782 nsteps = 2 

783 

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

785 force=False): 

786 

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

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

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

790 if self.store.config.stored_quantity == "rotation": 

791 self.gfmapping = [ 

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

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

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

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

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

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

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

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

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

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

802 ] 

803 else: 

804 self.gfmapping = [ 

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

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

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

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

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

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

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

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

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

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

815 ] 

816 

817 elif baseconf.qssp_version == 'ppeg2017': 

818 self.gfmapping = [ 

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

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

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

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

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

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

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

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

827 ] 

828 else: 

829 self.gfmapping = [ 

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

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

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

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

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

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

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

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

838 ] 

839 

840 if step == 0: 

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

842 else: 

843 if block_size is None: 

844 block_size = (1, 1, 51) 

845 

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

847 block_size = block_size[1:] 

848 

849 gf.builder.Builder.__init__( 

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

851 

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

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

854 conf.earthmodel_1d = self.store.config.earthmodel_1d 

855 deltat = self.store.config.deltat 

856 if 'time_window' not in shared: 

857 d = self.store.make_timing_params( 

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

859 force=force) 

860 

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

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

863 

864 shared['time_window'] = tmax - tmin 

865 shared['tstart'] = tmin 

866 

867 self.tstart = shared['tstart'] 

868 conf.time_window = shared['time_window'] 

869 

870 self.tmp = tmp 

871 if self.tmp is not None: 

872 util.ensuredir(self.tmp) 

873 

874 util.ensuredir(conf.gf_directory) 

875 

876 self.qssp_config = conf 

877 

878 def cleanup(self): 

879 self.store.close() 

880 

881 def work_block(self, iblock): 

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

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

884 self.get_block_extents(iblock) 

885 

886 rz = self.store.config.receiver_depth 

887 else: 

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

889 self.get_block_extents(iblock) 

890 

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

892 

893 conf = copy.deepcopy(self.qssp_config) 

894 

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

896 

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

898 logger.info( 

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

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

901 

902 return 

903 

904 logger.info( 

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

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

907 

908 tbeg = time.time() 

909 

910 runner = QSSPRunner(tmp=self.tmp) 

911 

912 conf.receiver_depth = rz/km 

913 conf.stored_quantity = self.store.config.stored_quantity 

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

915 dx = self.gf_config.distance_delta 

916 

917 if self.step == 0: 

918 distances = [firstx] 

919 else: 

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

921 

922 conf.receivers = [ 

923 QSSPReceiver( 

924 lat=90-d*cake.m2d, 

925 lon=180., 

926 tstart=self.tstart, 

927 distance=d) 

928 

929 for d in distances] 

930 

931 if self.step == 0: 

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

933 

934 gfs = [QSSPGreen( 

935 filename=gf_filename, 

936 depth=sz/km, 

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

938 

939 conf.greens_functions = gfs 

940 

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

942 

943 if self.step == 0: 

944 conf.sources = [QSSPSourceMT( 

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

946 lon=0.0, 

947 trise=trise, 

948 torigin=0.0)] 

949 

950 runner.run(conf) 

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

952 for s in glob.glob(gf_path): 

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

954 os.rename(s, d) 

955 

956 else: 

957 for mt, gfmap in self.gfmapping[ 

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

959 m = mt.m_up_south_east() 

960 

961 conf.sources = [QSSPSourceMT( 

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

963 lon=0.0, 

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

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

966 trise=trise, 

967 torigin=0.0)] 

968 

969 runner.run(conf) 

970 

971 rawtraces = runner.get_traces() 

972 

973 interrupted = [] 

974 

975 def signal_handler(signum, frame): 

976 interrupted.append(True) 

977 

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

979 self.store.lock() 

980 duplicate_inserts = 0 

981 try: 

982 for itr, tr in enumerate(rawtraces): 

983 if tr.channel in gfmap: 

984 

985 x = tr.meta['distance'] 

986 ig, factor = gfmap[tr.channel] 

987 

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

989 args = (sz, x, ig) 

990 else: 

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

992 

993 if conf.cut: 

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

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

996 if None in (tmin, tmax): 

997 continue 

998 

999 tr.chop(tmin, tmax) 

1000 

1001 tmin = tr.tmin 

1002 tmax = tr.tmax 

1003 

1004 if conf.fade: 

1005 ta, tb, tc, td = [ 

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

1007 for v in conf.fade] 

1008 

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

1010 continue 

1011 

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

1013 raise QSSPError( 

1014 'invalid fade configuration') 

1015 

1016 t = tr.get_xdata() 

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

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

1019 anti_fin = 1. - fin 

1020 anti_fout = 1. - fout 

1021 

1022 y = tr.ydata 

1023 

1024 sum_anti_fin = num.sum(anti_fin) 

1025 sum_anti_fout = num.sum(anti_fout) 

1026 

1027 if conf.nonzero_fade_in \ 

1028 and sum_anti_fin != 0.0: 

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

1030 else: 

1031 yin = 0.0 

1032 

1033 if conf.nonzero_fade_out \ 

1034 and sum_anti_fout != 0.0: 

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

1036 else: 

1037 yout = 0.0 

1038 

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

1040 

1041 if conf.relevel_with_fade_in: 

1042 y2 -= yin 

1043 

1044 tr.set_ydata(y2) 

1045 

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

1047 gf_tr.data *= factor 

1048 

1049 try: 

1050 self.store.put(args, gf_tr) 

1051 except gf.store.DuplicateInsert: 

1052 duplicate_inserts += 1 

1053 

1054 finally: 

1055 if duplicate_inserts: 

1056 logger.warn( 

1057 '%i insertions skipped (duplicates)' 

1058 % duplicate_inserts) 

1059 

1060 self.store.unlock() 

1061 signal.signal(signal.SIGINT, original) 

1062 

1063 if interrupted: 

1064 raise KeyboardInterrupt() 

1065 

1066 tend = time.time() 

1067 logger.info( 

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

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

1070 

1071 

1072km = 1000. 

1073 

1074 

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

1076 if variant is None: 

1077 variant = '2010beta' 

1078 

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

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

1081 

1082 qssp = QSSPConfig(qssp_version=variant) 

1083 if variant != 'ppeg2017': 

1084 qssp.time_region = ( 

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

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

1087 

1088 qssp.cut = ( 

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

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

1091 

1092 else: # variant == 'ppeg2017': 

1093 qssp.frequency_max = 0.5 

1094 qssp.time_region = [ 

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

1096 qssp.cut = [ 

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

1098 qssp.antialiasing_factor = 1.0e-10 

1099 qssp.toroidal_modes = False 

1100 qssp.cutoff_harmonic_degree_min = 2500 

1101 qssp.cutoff_harmonic_degree_max = 2500 

1102 qssp.crit_frequency_sge = 5.0 

1103 qssp.crit_harmonic_degree_sge = 50000 

1104 qssp.source_patch_radius = 10.0 

1105 qssp.bandpass_order = 6 

1106 qssp.bandpass_corner_low = 0.0 

1107 qssp.bandpass_corner_high = 0.125 

1108 

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

1110 if variant == 'ppeg2017': 

1111 quantity = 'acceleration' 

1112 else: 

1113 quantity = None 

1114 

1115 if variant == 'ppeg2017': 

1116 sample_rate = 4.0 

1117 else: 

1118 sample_rate = 0.2 

1119 

1120 d = dict( 

1121 id=store_id, 

1122 ncomponents=10, 

1123 component_scheme='elastic10', 

1124 stored_quantity=quantity, 

1125 sample_rate=sample_rate, 

1126 receiver_depth=0*km, 

1127 source_depth_min=10*km, 

1128 source_depth_max=20*km, 

1129 source_depth_delta=10*km, 

1130 distance_min=100*km, 

1131 distance_max=1000*km, 

1132 distance_delta=10*km, 

1133 earthmodel_1d=cake.load_model(), 

1134 modelling_code_id='qssp', 

1135 tabulated_phases=[ 

1136 gf.meta.TPDef( 

1137 id='begin', 

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

1139 gf.meta.TPDef( 

1140 id='end', 

1141 definition='2.5'), 

1142 gf.meta.TPDef( 

1143 id='P', 

1144 definition='!P'), 

1145 gf.meta.TPDef( 

1146 id='S', 

1147 definition='!S'), 

1148 gf.meta.TPDef( 

1149 id='p', 

1150 definition='!p'), 

1151 gf.meta.TPDef( 

1152 id='s', 

1153 definition='!s')]) 

1154 

1155 if config_params is not None: 

1156 d.update(config_params) 

1157 

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

1159 

1160 config.validate() 

1161 return gf.store.Store.create_editables( 

1162 store_dir, 

1163 config=config, 

1164 extra={'qssp': qssp}) 

1165 

1166 

1167def build( 

1168 store_dir, 

1169 force=False, 

1170 nworkers=None, 

1171 continue_=False, 

1172 step=None, 

1173 iblock=None): 

1174 

1175 return QSSPGFBuilder.build( 

1176 store_dir, force=force, nworkers=nworkers, continue_=continue_, 

1177 step=step, iblock=iblock) 

1178 

1179 

1180def get_conf( 

1181 store_dir, 

1182 force=False, 

1183 nworkers=None, 

1184 continue_=False, 

1185 step=None, 

1186 iblock=None): 

1187 

1188 return QSSPGFBuilder.get_conf()