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

450 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2024-03-07 11:54 +0000

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_displacement': 

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_displacement': 

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.component_scheme == 'rotational8': 

790 self.gfmapping = [ 

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

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

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

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

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

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

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

798 {'_rota_e': (4, -1)}), 

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

800 {'_rota_e': (7, -1)}), 

801 ] 

802 elif self.store.config.component_scheme == 'elastic10': 

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 elif self.store.config.component_scheme == 'elastic8': 

816 self.gfmapping = [ 

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

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

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

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

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

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

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

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

825 ] 

826 

827 elif baseconf.qssp_version == 'ppeg2017': 

828 self.gfmapping = [ 

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

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

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

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

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

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

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

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

837 ] 

838 else: 

839 self.gfmapping = [ 

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

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

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

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

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

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

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

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

848 ] 

849 

850 if step == 0: 

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

852 else: 

853 if block_size is None: 

854 block_size = (1, 1, 51) 

855 

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

857 block_size = block_size[1:] 

858 

859 gf.builder.Builder.__init__( 

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

861 

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

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

864 conf.earthmodel_1d = self.store.config.earthmodel_1d 

865 deltat = self.store.config.deltat 

866 if 'time_window' not in shared: 

867 d = self.store.make_timing_params( 

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

869 force=force) 

870 

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

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

873 

874 shared['time_window'] = tmax - tmin 

875 shared['tstart'] = tmin 

876 

877 self.tstart = shared['tstart'] 

878 conf.time_window = shared['time_window'] 

879 

880 self.tmp = tmp 

881 if self.tmp is not None: 

882 util.ensuredir(self.tmp) 

883 

884 util.ensuredir(conf.gf_directory) 

885 

886 self.qssp_config = conf 

887 

888 def cleanup(self): 

889 self.store.close() 

890 

891 def work_block(self, iblock): 

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

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

894 self.get_block_extents(iblock) 

895 

896 rz = self.store.config.receiver_depth 

897 else: 

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

899 self.get_block_extents(iblock) 

900 

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

902 

903 conf = copy.deepcopy(self.qssp_config) 

904 

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

906 

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

908 logger.info( 

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

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

911 

912 return 

913 

914 logger.info( 

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

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

917 

918 tbeg = time.time() 

919 

920 runner = QSSPRunner(tmp=self.tmp) 

921 

922 conf.receiver_depth = rz/km 

923 

924 component_scheme = self.store.config.component_scheme 

925 stored_quantity = self.store.config.effective_stored_quantity 

926 

927 if (component_scheme, stored_quantity) not in [ 

928 ('rotational8', 'rotation_displacement'), 

929 ('elastic10', 'displacement'), 

930 ('elastic8', 'displacement')]: 

931 

932 raise QSSPError( 

933 'componentent_scheme "%s" is inconsistent with ' 

934 'stored_quantity "%s"' % (component_scheme, stored_quantity)) 

935 

936 conf.stored_quantity = stored_quantity 

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

938 dx = self.gf_config.distance_delta 

939 

940 if self.step == 0: 

941 distances = [firstx] 

942 else: 

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

944 

945 conf.receivers = [ 

946 QSSPReceiver( 

947 lat=90-d*cake.m2d, 

948 lon=180., 

949 tstart=self.tstart, 

950 distance=d) 

951 

952 for d in distances] 

953 

954 if self.step == 0: 

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

956 

957 gfs = [QSSPGreen( 

958 filename=gf_filename, 

959 depth=sz/km, 

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

961 

962 conf.greens_functions = gfs 

963 

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

965 

966 if self.step == 0: 

967 conf.sources = [QSSPSourceMT( 

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

969 lon=0.0, 

970 trise=trise, 

971 torigin=0.0)] 

972 

973 runner.run(conf) 

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

975 for s in glob.glob(gf_path): 

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

977 os.rename(s, d) 

978 

979 else: 

980 for mt, gfmap in self.gfmapping: 

981 m = mt.m_up_south_east() 

982 

983 conf.sources = [QSSPSourceMT( 

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

985 lon=0.0, 

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

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

988 trise=trise, 

989 torigin=0.0)] 

990 

991 runner.run(conf) 

992 

993 rawtraces = runner.get_traces() 

994 

995 interrupted = [] 

996 

997 def signal_handler(signum, frame): 

998 interrupted.append(True) 

999 

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

1001 self.store.lock() 

1002 duplicate_inserts = 0 

1003 try: 

1004 for itr, tr in enumerate(rawtraces): 

1005 if tr.channel in gfmap: 

1006 

1007 x = tr.meta['distance'] 

1008 ig, factor = gfmap[tr.channel] 

1009 

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

1011 args = (sz, x, ig) 

1012 else: 

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

1014 

1015 if conf.cut: 

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

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

1018 if None in (tmin, tmax): 

1019 continue 

1020 

1021 tr.chop(tmin, tmax) 

1022 

1023 tmin = tr.tmin 

1024 tmax = tr.tmax 

1025 

1026 if conf.fade: 

1027 ta, tb, tc, td = [ 

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

1029 for v in conf.fade] 

1030 

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

1032 continue 

1033 

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

1035 raise QSSPError( 

1036 'invalid fade configuration') 

1037 

1038 t = tr.get_xdata() 

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

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

1041 anti_fin = 1. - fin 

1042 anti_fout = 1. - fout 

1043 

1044 y = tr.ydata 

1045 

1046 sum_anti_fin = num.sum(anti_fin) 

1047 sum_anti_fout = num.sum(anti_fout) 

1048 

1049 if conf.nonzero_fade_in \ 

1050 and sum_anti_fin != 0.0: 

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

1052 else: 

1053 yin = 0.0 

1054 

1055 if conf.nonzero_fade_out \ 

1056 and sum_anti_fout != 0.0: 

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

1058 else: 

1059 yout = 0.0 

1060 

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

1062 

1063 if conf.relevel_with_fade_in: 

1064 y2 -= yin 

1065 

1066 tr.set_ydata(y2) 

1067 

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

1069 gf_tr.data *= factor 

1070 

1071 try: 

1072 self.store.put(args, gf_tr) 

1073 except gf.store.DuplicateInsert: 

1074 duplicate_inserts += 1 

1075 

1076 finally: 

1077 if duplicate_inserts: 

1078 logger.warning( 

1079 '%i insertions skipped (duplicates)' 

1080 % duplicate_inserts) 

1081 

1082 self.store.unlock() 

1083 signal.signal(signal.SIGINT, original) 

1084 

1085 if interrupted: 

1086 raise KeyboardInterrupt() 

1087 

1088 tend = time.time() 

1089 logger.info( 

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

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

1092 

1093 

1094km = 1000. 

1095 

1096 

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

1098 if variant is None: 

1099 variant = '2010beta' 

1100 

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

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

1103 

1104 qssp = QSSPConfig(qssp_version=variant) 

1105 if variant != 'ppeg2017': 

1106 qssp.time_region = ( 

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

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

1109 

1110 qssp.cut = ( 

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

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

1113 

1114 else: # variant == 'ppeg2017': 

1115 qssp.frequency_max = 0.5 

1116 qssp.time_region = [ 

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

1118 qssp.cut = [ 

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

1120 qssp.antialiasing_factor = 1.0e-10 

1121 qssp.toroidal_modes = False 

1122 qssp.cutoff_harmonic_degree_min = 2500 

1123 qssp.cutoff_harmonic_degree_max = 2500 

1124 qssp.crit_frequency_sge = 5.0 

1125 qssp.crit_harmonic_degree_sge = 50000 

1126 qssp.source_patch_radius = 10.0 

1127 qssp.bandpass_order = 6 

1128 qssp.bandpass_corner_low = 0.0 

1129 qssp.bandpass_corner_high = 0.125 

1130 

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

1132 if variant == 'ppeg2017': 

1133 quantity = 'acceleration' 

1134 else: 

1135 quantity = None 

1136 

1137 if variant == 'ppeg2017': 

1138 sample_rate = 4.0 

1139 else: 

1140 sample_rate = 0.2 

1141 

1142 d = dict( 

1143 id=store_id, 

1144 ncomponents=10, 

1145 component_scheme='elastic10', 

1146 stored_quantity=quantity, 

1147 sample_rate=sample_rate, 

1148 receiver_depth=0*km, 

1149 source_depth_min=10*km, 

1150 source_depth_max=20*km, 

1151 source_depth_delta=10*km, 

1152 distance_min=100*km, 

1153 distance_max=1000*km, 

1154 distance_delta=10*km, 

1155 earthmodel_1d=cake.load_model(), 

1156 modelling_code_id='qssp', 

1157 tabulated_phases=[ 

1158 gf.meta.TPDef( 

1159 id='begin', 

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

1161 gf.meta.TPDef( 

1162 id='end', 

1163 definition='2.5'), 

1164 gf.meta.TPDef( 

1165 id='P', 

1166 definition='!P'), 

1167 gf.meta.TPDef( 

1168 id='S', 

1169 definition='!S'), 

1170 gf.meta.TPDef( 

1171 id='p', 

1172 definition='!p'), 

1173 gf.meta.TPDef( 

1174 id='s', 

1175 definition='!s')]) 

1176 

1177 if config_params is not None: 

1178 d.update(config_params) 

1179 

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

1181 

1182 config.validate() 

1183 return gf.store.Store.create_editables( 

1184 store_dir, 

1185 config=config, 

1186 extra={'qssp': qssp}) 

1187 

1188 

1189def build( 

1190 store_dir, 

1191 force=False, 

1192 nworkers=None, 

1193 continue_=False, 

1194 step=None, 

1195 iblock=None): 

1196 

1197 return QSSPGFBuilder.build( 

1198 store_dir, force=force, nworkers=nworkers, continue_=continue_, 

1199 step=step, iblock=iblock) 

1200 

1201 

1202def get_conf( 

1203 store_dir, 

1204 force=False, 

1205 nworkers=None, 

1206 continue_=False, 

1207 step=None, 

1208 iblock=None): 

1209 

1210 return QSSPGFBuilder.get_conf()