Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/fomosto/report/report_main.py: 66%

1002 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 os 

7import sys 

8import shutil 

9import subprocess 

10import re 

11import io 

12import base64 

13import datetime 

14from collections import OrderedDict 

15from tempfile import NamedTemporaryFile, mkdtemp 

16from string import Template 

17 

18import numpy as np 

19from matplotlib import pyplot as plt 

20from matplotlib import cm, transforms 

21 

22from pyrocko import gf, trace, cake, util, plot 

23from pyrocko.response import DifferentiationResponse 

24from pyrocko.plot import beachball 

25from pyrocko.guts import load, Object, String, List, Float, Int, Bool, Dict 

26from pyrocko.gf import Source, Target 

27from pyrocko.gf.meta import OutOfBounds 

28 

29from jinja2 import Environment, PackageLoader 

30 

31guts_prefix = 'gft' 

32ex_path = os.path.dirname(os.path.abspath(sys.argv[0])) 

33ja_latex_env = Environment(block_start_string='\\BLOCK{', 

34 block_end_string='}', 

35 variable_start_string='\\VAR{', 

36 variable_end_string='}', 

37 comment_start_string='\\%{', 

38 comment_end_string='}', 

39 line_statement_prefix='\\-', 

40 line_comment_prefix='%%', 

41 trim_blocks=True, 

42 loader=PackageLoader('pyrocko', 

43 'data/fomosto_report')) 

44 

45 

46class FomostoReportError(Exception): 

47 pass 

48 

49 

50class GreensFunctionError(FomostoReportError): 

51 pass 

52 

53 

54class FilterFrequencyError(GreensFunctionError): 

55 

56 def __init__(self, frequency): 

57 Exception.__init__(self) 

58 self.frequency = frequency 

59 

60 def __str__(self): 

61 return 'Cannot set {0} frequency to both an absolute and' \ 

62 ' revlative value.'.format(self.frequency) 

63 

64 

65class SourceError(GreensFunctionError): 

66 

67 def __init__(self, typ): 

68 Exception.__init__(self) 

69 self.type = typ 

70 

71 def __str__(self): 

72 return 'Source type not currently supported: {0}'.format(self.type) 

73 

74 

75class SensorArray(Target): 

76 

77 distance_min = Float.T() 

78 distance_max = Float.T() 

79 strike = Float.T() 

80 sensor_count = Int.T(default=50) 

81 

82 __target_allowed = dir(Target) 

83 __sensorarray_allowed = ['distance_min', 'distance_max', 'strike', 

84 'sensor_count', 'change_dists'] 

85 

86 def __init__(self, **kwargs): 

87 Object.__init__(self, **kwargs) 

88 

89 self.validate_args(kwargs) 

90 self.__build(**kwargs) 

91 

92 @classmethod 

93 def validate_args(cls, kwargs): 

94 ks = [] 

95 for k in kwargs: 

96 if not (k in cls.__target_allowed or 

97 k in cls.__sensorarray_allowed): 

98 ks.append(k) 

99 

100 for i in ks: 

101 del kwargs[i] 

102 

103 @classmethod 

104 def __validate_args_target(cls, kwargs): 

105 ks = [] 

106 for k in kwargs: 

107 if k not in cls.__target_allowed: 

108 ks.append(k) 

109 

110 for i in ks: 

111 del kwargs[i] 

112 

113 def __build(self, **kwargs): 

114 if 'quantity' not in kwargs: 

115 kwargs['quantity'] = 'displacement' 

116 dists = np.linspace(self.distance_min, self.distance_max, 

117 self.sensor_count) 

118 self.__validate_args_target(kwargs) 

119 self.sensors = [ 

120 Target(north_shift=float(np.cos(np.radians(self.strike)) * dist), 

121 east_shift=float(np.sin(np.radians(self.strike)) * dist), 

122 **kwargs) 

123 for dist in dists] 

124 

125 

126class GreensFunctionTest(Object): 

127 

128 __valid_properties = ['store_dir', 'pdf_dir', 'output_format', 

129 'lowpass_frequency', 'rel_lowpass_frequency', 

130 'highpass_frequency', 'rel_highpass_frequency', 

131 'filter_order', 

132 'plot_velocity', 'plot_everything'] 

133 __notesize = 7.45 

134 __scalelist = [1, 5, 9.5, 19, 29] 

135 __has_phase = True 

136 

137 store_dir = String.T() 

138 pdf_dir = String.T(default=os.path.expanduser('~')) 

139 output_format = String.T(optional=True, default='pdf') 

140 lowpass_frequency = Float.T(optional=True) 

141 highpass_frequency = Float.T(optional=True) 

142 rel_lowpass_frequency = Float.T(optional=True) 

143 rel_highpass_frequency = Float.T(optional=True) 

144 filter_order = Int.T(default=4) 

145 sensor_count = Int.T(default=50) 

146 plot_velocity = Bool.T(default=False) 

147 plot_everything = Bool.T(default=True) 

148 src_ids = List.T(String.T()) 

149 sen_ids = List.T(String.T()) 

150 sources = Dict.T(String.T(), Source.T()) 

151 sensors = Dict.T(String.T(), SensorArray.T()) 

152 trace_configs = List.T(List.T(String.T()), optional=True) 

153 

154 @classmethod 

155 def __get_valid_arguments(cls, args): 

156 tdict = {} 

157 for k in args: 

158 if k in cls.__valid_properties: 

159 tdict[k] = args[k] 

160 return tdict 

161 

162 def __init__(self, store_dir, **kwargs): 

163 Object.__init__(self, store_dir=store_dir, **kwargs) 

164 

165 if self.lowpass_frequency and self.rel_lowpass_frequency: 

166 raise FilterFrequencyError('lowpass') 

167 if self.highpass_frequency and self.rel_highpass_frequency: 

168 raise FilterFrequencyError('highpass') 

169 

170 if self.store_dir[-1] != '/': 

171 self.store_dir += '/' 

172 self.engine = gf.LocalEngine(store_dirs=[self.store_dir]) 

173 ids = self.engine.get_store_ids() 

174 self.store_id = ids[0] 

175 self.store = self.engine.get_store(self.store_id) 

176 

177 if self.rel_lowpass_frequency is not None: 

178 self.lowpass_frequency = self.store.config.sample_rate * \ 

179 self.rel_lowpass_frequency 

180 

181 if self.rel_highpass_frequency is not None: 

182 self.highpass_frequency = self.store.config.sample_rate * \ 

183 self.rel_highpass_frequency 

184 

185 self.rel_lowpass_frequency = None 

186 self.rel_highpass_frequency = None 

187 

188 self.phases = '|'.join([x.id 

189 for x in self.store.config.tabulated_phases]) 

190 if self.phases == '': 

191 self.phase_ratio_string = r'\color{red}' \ 

192 '{warning: store has no tabulated phases}' 

193 elif 'begin' in self.phases: 

194 self.phase_ratio_string = 'begin' 

195 else: 

196 self.phase_ratio_string = 'first({0})'.format(self.phases) 

197 

198 if len(self.src_ids) == 0 and len(self.sources) > 0: 

199 self.src_ids = sorted(self.sources) 

200 if len(self.sen_ids) == 0 and len(self.sensors) > 0: 

201 self.sen_ids = sorted(self.sensors) 

202 

203 self.traces = OrderedDict() 

204 

205 if self.pdf_dir is None: 

206 self.pdf_dir = os.path.expanduser('~') 

207 if self.pdf_dir[-1] != '/': 

208 self.pdf_dir += '/' 

209 self.pdf_name = '{0}_{1}'.format( 

210 self.store_id, self.getFrequencyString(None)) 

211 util.setup_logging() 

212 

213 self.temp_dir = mkdtemp(prefix='gft_') 

214 self.message = None 

215 self.changed_depth = False 

216 self.changed_dist_min = False 

217 self.changed_dist_max = False 

218 

219 def __addToMessage(self, msg): 

220 if self.message is None: 

221 self.message = '\nmessage(s) for store {0}:\n'. \ 

222 format(self.store_id) 

223 self.message += msg + '\n' 

224 

225 def __printMessage(self): 

226 if self.message is not None: 

227 print(self.message) 

228 

229 def cleanup(self): 

230 shutil.rmtree(self.temp_dir) 

231 for i in ['.aux', '.log', '.out', '.toc']: 

232 path = '{0}{1}{2}'.format(self.pdf_dir, self.pdf_name, i) 

233 if os.path.exists(path): 

234 os.remove(path) 

235 

236 def getSourceString(self, sid): 

237 if sid not in self.sources: 

238 return '' 

239 src = self.sources[sid] 

240 if isinstance(src, gf.DCSource): 

241 typ = 'Double Couple' 

242 sstr = 'Source Type: {0}, Strike: {1:g}, Dip: {2:g}, Rake: {3:g}' \ 

243 .format(typ, src.strike, src.dip, src.rake) 

244 elif isinstance(src, gf.RectangularExplosionSource): 

245 typ = 'Explosion' 

246 sstr = 'Source Type: {0}, Strike: {1:g}, Dip: {2:g}'.format( 

247 typ, src.strike, src.dip) 

248 else: 

249 typ = '{0}'.format(type(src)).split('.')[-1].split("'")[0] 

250 sstr = 'Source Type: {0}, see config page'.format(typ) 

251 return sstr 

252 

253 def getSensorString(self, sid): 

254 if sid not in self.sensors: 

255 return '' 

256 sen_ar = self.sensors[sid] 

257 sensors = sen_ar.sensors 

258 

259 targ = sensors[0] 

260 senstr = 'Sensor Strike: {0:g}, Azimuth: {1:g}, Dip: {2:g}'. \ 

261 format(sen_ar.strike, targ.azimuth, targ.dip) 

262 return senstr 

263 

264 def getFrequencyString(self, tid): 

265 lpf = self.lowpass_frequency 

266 hpf = self.highpass_frequency 

267 if tid is None: 

268 if lpf and hpf: 

269 return '{0:.4g}-{1:.4g}Hz'.format(lpf, hpf) 

270 elif lpf: 

271 return 'lowpass {0:.4g}Hz'.format(lpf) 

272 elif hpf: 

273 return 'highpass {0:.4g}Hz'.format(hpf) 

274 else: 

275 return 'unfiltered' 

276 

277 if tid not in self.traces: 

278 return '' 

279 

280 tdict = self.traces[tid] 

281 lpfa = tdict['lowpass_applied'] 

282 hpfa = tdict['highpass_applied'] 

283 if lpf and hpf and lpfa and hpfa: 

284 return '{0:.4g}-{1:.4g} [Hz]'.format(lpf, hpf) 

285 elif lpf and lpfa: 

286 return 'lowpass frequency: {0:.4g} [Hz]'.format(lpf) 

287 elif hpf and hpfa: 

288 return 'highpass frequency: {0:.4g} [Hz]'.format(hpf) 

289 else: 

290 return 'unfiltered' 

291 

292 def getAnnotateString(self, src_id, sen_id, fac=None): 

293 tstr = "Green's Function: {0}\n{1}\n{2}".format( 

294 self.store_id, self.getSourceString(src_id), 

295 self.getSensorString(sen_id)) 

296 if fac is not None: 

297 tstr += '\nAmplitude Scale Factor: {0:0.4g} [microns/km]'. \ 

298 format(fac * 1e-9 / 1e3) 

299 return tstr 

300 

301 def createSource(self, typ, depth, strike, dip, rake=None, lat=0., lon=0., 

302 north_shift=0., east_shift=0., change_dists=True, 

303 **kwargs): 

304 stCfg = self.store.config 

305 stMin = stCfg.source_depth_min 

306 stMax = stCfg.source_depth_max 

307 if depth is None: 

308 depth = (stMax + stMin) / 2. 

309 elif depth < stMin or depth > stMax: 

310 if not change_dists: 

311 raise OutOfBounds('Source depth.') 

312 diff = abs(stMin - depth) 

313 if diff < stMax - depth: 

314 depth = stMin 

315 else: 

316 depth = stMax 

317 if not self.changed_depth: 

318 self.changed_depth = True 

319 self.__addToMessage( 

320 'Source depth out of bounds. Changed to: {0:g}m.'. 

321 format(depth)) 

322 

323 if typ.upper() in ('EX' 'EXPLOSION'): 

324 src = gf.RectangularExplosionSource( 

325 lat=lat, 

326 lon=lon, 

327 north_shift=north_shift, 

328 east_shift=east_shift, 

329 depth=depth, 

330 strike=strike, 

331 dip=dip) 

332 elif typ.upper() in ('DC' 'DOUBLE COUPLE'): 

333 src = gf.DCSource( 

334 lat=lat, 

335 lon=lon, 

336 north_shift=north_shift, 

337 east_shift=east_shift, 

338 depth=depth, 

339 strike=strike, 

340 dip=dip, 

341 rake=rake) 

342 else: 

343 raise SourceError(typ) 

344 

345 srcstr = '{0}.{1:.2g}.{2:.2g}.{3:.2g}'.format( 

346 typ[0:2], depth, strike, dip) 

347 tstr = srcstr 

348 tint = 96 

349 while tstr in self.sources: 

350 tint += 1 

351 tstr = srcstr + chr(tint) 

352 

353 self.sources[tstr] = src 

354 self.src_ids.append(tstr) 

355 return tstr 

356 

357 def createSensors(self, distance_min=None, distance_max=None, 

358 change_dists=True, **kwargs): 

359 stCfg = self.store.config 

360 stMin = stCfg.distance_min 

361 stMax = stCfg.distance_max 

362 if distance_min is None: 

363 distance_min = stMin 

364 elif distance_min < stMin or distance_min > stMax: 

365 if not change_dists: 

366 raise OutOfBounds('Minimum sensor distance.') 

367 distance_min = stMin 

368 if not self.changed_dist_min: 

369 self.changed_dist_min = True 

370 self.__addToMessage( 

371 'Sensor minimum distance out of bounds. Changed to allowed' 

372 ' minimum {0:g}m.'.format(stMin)) 

373 

374 if distance_max is None: 

375 distance_max = stMax 

376 elif (distance_max > stMax or distance_max < stMin): 

377 if not change_dists: 

378 raise OutOfBounds('Maximum sensor distance') 

379 distance_max = stMax 

380 if not self.changed_dist_max: 

381 self.changed_dist_max = True 

382 self.__addToMessage( 

383 'Sensor maximum distance out of bounds. Changed to allowed' 

384 ' maximum {0:g}m.'.format(stMax)) 

385 

386 if change_dists and distance_min == distance_max: 

387 distance_min = stMin 

388 distance_max = stMax 

389 self.__addToMessage( 

390 'Sensor minimum and maximum distance equal. Changed to' 

391 ' allowed minimum {0:g}m and maximum {1:g}m.'.format(stMin, 

392 stMax)) 

393 

394 sen_ar = SensorArray(distance_min=distance_min, 

395 distance_max=distance_max, **kwargs) 

396 senstr = '{0}.{1:.2g}'.format('.'.join(sen_ar.codes), sen_ar.strike) 

397 tstr = senstr 

398 tint = 96 

399 while tstr in self.sensors: 

400 tint += 1 

401 tstr = senstr + chr(tint) 

402 

403 self.sensors[tstr] = sen_ar 

404 self.sen_ids.append(tstr) 

405 return tstr 

406 

407 def createDisplacementTraces(self, src_id='all', sen_id='all'): 

408 for sr_id in self.src_ids: 

409 if sr_id not in self.sources: 

410 continue 

411 if not (src_id == 'all' or sr_id == src_id): 

412 continue 

413 for sn_id in self.sen_ids: 

414 if sn_id not in self.sensors: 

415 continue 

416 if not (sen_id == 'all' or sn_id == sen_id): 

417 continue 

418 try: 

419 response = self.engine.process( 

420 self.sources[sr_id], 

421 self.sensors[sn_id].sensors) 

422 trcs = response.pyrocko_traces() 

423 tstr = '{0}|{1}'.format(sr_id, sn_id) 

424 if tstr not in self.traces: 

425 self.traces[tstr] = {} 

426 tdict = self.traces[tstr] 

427 tdict['displacement_traces'] = trcs 

428 mina, maxa, minb, maxb, ratio = \ 

429 self.__tracesMinMax(trcs, sr_id, sn_id) 

430 if ratio != 0.: 

431 tdict['displacement_spectra'] = [ 

432 trc.spectrum() for trc in trcs] 

433 tdict['lowpass_applied'] = False 

434 tdict['highpass_applied'] = False 

435 tdict['displacement_ratio'] = ratio 

436 tdict['displacement_scale'] = max(abs(mina), abs(maxa)) 

437 except IndexError: 

438 self.__addToMessage( 

439 'warning: IndexError: no traces created for' 

440 ' source-sensor combination: {0} - {1}. Try increasing' 

441 ' the sensor minimum distance.'.format( 

442 sr_id, sn_id)) 

443 

444 def createVelocityTraces(self, trc_id='all'): 

445 for tid in self.traces: 

446 ids = trc_id.split('|') 

447 if len(ids) == 2: 

448 src_id, sen_id = ids 

449 else: 

450 src_id = sen_id = None 

451 if not (trc_id == 'all' or 

452 ((src_id == 'all' or src_id in self.sources) and 

453 (sen_id == 'all' or sen_id in self.sensors))): 

454 continue 

455 tdict = self.traces[tid] 

456 if 'displacement_traces' not in tdict: 

457 continue 

458 trcs = [trc.transfer( 

459 transfer_function=DifferentiationResponse()) 

460 for trc in tdict['displacement_traces']] 

461 tdict['velocity_traces'] = trcs 

462 src_id, sen_id = tid.split('|') 

463 mina, maxa, minb, maxb, ratio = \ 

464 self.__tracesMinMax(trcs, src_id, sen_id) 

465 if ratio != 0.: 

466 tdict['velocity_spectra'] = [ 

467 trc.spectrum() for trc in trcs] 

468 tdict['velocity_ratio'] = ratio 

469 tdict['velocity_scale'] = max(abs(mina), abs(maxa)) 

470 

471 def getPhaseArrivals(self, trc_id='all', phase_ids='all'): 

472 for tid in self.traces: 

473 if trc_id == 'all' or trc_id == tid: 

474 self.__setPhaseArrivals(tid, phase_ids) 

475 

476 def __setPhaseArrivals(self, trc_id, phase_ids='all'): 

477 st = self.store 

478 src_id, sen_id = trc_id.split('|') 

479 src = self.sources[src_id] 

480 tdict = self.traces[trc_id] 

481 if 'phase_arrivals' not in tdict: 

482 tdict['phase_arrivals'] = {} 

483 phdict = tdict['phase_arrivals'] 

484 for pid in self.phases.split('|'): 

485 if pid == '': 

486 continue 

487 if phase_ids == 'all' or pid == phase_ids: 

488 dists = [] 

489 times = [] 

490 for sen in self.sensors[sen_id].sensors: 

491 time = st.t(pid, src, sen) 

492 if time is None: 

493 continue 

494 dists.append(src.distance_to(sen)) 

495 times.append(time) 

496 

497 if len(times) > 0: 

498 phdict[pid] = (times, dists) 

499 

500 def applyFrequencyFilters(self, trc_id='all'): 

501 dispname = 'displacement_traces' 

502 velname = 'velocity_traces' 

503 for tid in self.traces: 

504 if trc_id == 'all' or trc_id == tid: 

505 tdict = self.traces[tid] 

506 if self.lowpass_frequency: 

507 if dispname in tdict: 

508 for trc in tdict[dispname]: 

509 trc.lowpass(self.filter_order, 

510 self.lowpass_frequency, demean=False) 

511 tdict['lowpass_applied'] = True 

512 

513 if velname in tdict: 

514 for trc in tdict[velname]: 

515 trc.lowpass(self.filter_order, 

516 self.lowpass_frequency, demean=False) 

517 tdict['lowpass_applied'] = True 

518 

519 if self.highpass_frequency: 

520 if dispname in tdict: 

521 for trc in tdict[dispname]: 

522 trc.highpass(self.filter_order, 

523 self.highpass_frequency, demean=False) 

524 tdict['highpass_applied'] = True 

525 

526 if velname in tdict: 

527 for trc in tdict[velname]: 

528 trc.highpass(self.filter_order, 

529 self.highpass_frequency, demean=False) 

530 tdict['highpass_applied'] = True 

531 

532 sr_id, sn_id = tid.split('|') 

533 if dispname in tdict: 

534 trcs = tdict[dispname] 

535 mina, maxa, minb, maxb, ratio = \ 

536 self.__tracesMinMax(trcs, sr_id, sn_id) 

537 tdict['displacement_ratio'] = ratio 

538 tdict['displacement_scale'] = max(abs(mina), abs(maxa)) 

539 if velname in tdict: 

540 trcs = tdict[velname] 

541 mina, maxa, minb, maxb, ratio = \ 

542 self.__tracesMinMax(trcs, sr_id, sn_id) 

543 tdict['velocity_ratio'] = ratio 

544 tdict['velocity_scale'] = max(abs(mina), abs(maxa)) 

545 

546 def __createOutputDoc(self, artefacts, chapters, 

547 gft2=None, together=False): 

548 

549 str_id = self.store_id 

550 is_tex = self.output_format == 'pdf' 

551 if is_tex: 

552 file_type = 'tex' 

553 dir = self.temp_dir 

554 fstr_id = self.__formatLatexString(str_id) 

555 else: 

556 file_type = self.output_format 

557 dir = self.pdf_dir 

558 fstr_id = str_id 

559 

560 temp = ja_latex_env.get_template('gfreport.{0}'.format(file_type)) 

561 out = '{0}/{1}.{2}'.format(dir, self.pdf_name, file_type) 

562 config = self.store.config.dump() 

563 

564 if together: 

565 tpath = self.__createVelocityProfile(gft2) 

566 if is_tex: 

567 img_ttl = r'{0} \& {1}'.format( 

568 fstr_id, gft2.__formatLatexString(gft2.store_id)) 

569 else: 

570 img_ttl = r'{0} & {1}'.format(str_id, gft2.store_id) 

571 else: 

572 tpath = self.__createVelocityProfile() 

573 img_ttl = fstr_id 

574 

575 info = [(fstr_id, config, tpath, img_ttl)] 

576 

577 if gft2 is not None: 

578 if is_tex: 

579 fstr_id = self.__formatLatexString(gft2.store_id) 

580 str_id = r'{0} \& {1}'.format(str_id, gft2.store_id) 

581 else: 

582 fstr_id = gft2.store_id 

583 str_id = r'{0} & {1}'.format(str_id, gft2.store_id) 

584 

585 config = gft2.store.config.dump() 

586 if together: 

587 tpath = '' 

588 else: 

589 tpath = gft2.__createVelocityProfile() 

590 

591 info += [(fstr_id, config, tpath, fstr_id)] 

592 

593 with open(out, 'w') as f: 

594 if is_tex: 

595 f.write(temp.render( 

596 rpt_id=self.__formatLatexString(str_id), str_info=info, 

597 artefacts=artefacts, chapters=chapters, 

598 config=config, headings='headings')) 

599 elif file_type == 'html': 

600 f.write(temp.render( 

601 rpt_id=str_id, str_info=info, artefacts=artefacts, 

602 chapters=chapters, config=config, 

603 date=datetime.datetime.now().strftime('%B %d, %Y'))) 

604 

605 if is_tex: 

606 pro_call = ['pdflatex', '-output-directory', self.pdf_dir, out, 

607 '-interaction', 'nonstop'] 

608 try: 

609 subprocess.call(pro_call) 

610 subprocess.call(pro_call) 

611 subprocess.call(pro_call) 

612 except OSError: 

613 raise FomostoReportError( 

614 'Cannot run "pdflatex" executable. Is it installed?') 

615 

616 self.cleanup() 

617 if gft2 is not None: 

618 gft2.cleanup() 

619 

620 @staticmethod 

621 def __formatLatexString(string): 

622 rep = {'_': r'\_', r'\n': r'\\', '|': r'\textbar '} 

623 rep = {re.escape(k): v for k, v in rep.items()} 

624 pat = re.compile('|'.join(rep.keys())) 

625 return pat.sub(lambda m: rep[re.escape(m.group(0))], string) 

626 

627 def createOutputDoc(self, *trc_ids): 

628 trcs = self.traces 

629 if len(trc_ids) == 0: 

630 trc_ids = trcs.keys() 

631 artefacts = self.__getArtefactPageInfo(trc_ids) 

632 chapters = [] 

633 if self.plot_everything: 

634 for tid in trc_ids: 

635 if tid not in trcs: 

636 continue 

637 src_id, sen_id = tid.split('|') 

638 ttl = '{0}, {1}'.format( 

639 self.getSourceString(src_id), 

640 self.getSensorString(sen_id)) 

641 tdict = trcs[tid] 

642 img_data = [] 

643 href = r'\hypertarget{${trc_id}|${str_id}|${type}}{}' 

644 trcname = 'displacement_traces' 

645 if trcname in tdict: 

646 hstr = Template(href).substitute( 

647 trc_id=tid, str_id=self.store_id, type='d') 

648 figs = self.__createTraceFigures(tid, trcname) 

649 fig = figs[0] 

650 img_data.extend([(self.__getFigureTitle(fig), hstr, 

651 self.__saveTempFigure(fig))]) 

652 img_data.extend([('', '', self.__saveTempFigure(x)) 

653 for x in figs[1:]]) 

654 

655 if 'displacement_spectra' in tdict: 

656 fig = self.__createMaxAmpFigure(tid, trcname) 

657 img_data.append((self.__getFigureTitle(fig), '', 

658 self.__saveTempFigure(fig))) 

659 

660 fig = self.__createSpectraFigure(tid, 

661 'displacement_spectra') 

662 img_data.append((self.__getFigureTitle(fig), '', 

663 self.__saveTempFigure(fig))) 

664 

665 trcname = 'velocity_traces' 

666 if self.plot_velocity and trcname in tdict: 

667 hstr = Template(href).substitute( 

668 trc_id=tid, str_id=self.store_id, type='v') 

669 figs = self.__createTraceFigures(tid, trcname) 

670 fig = figs[0] 

671 img_data.extend([(self.__getFigureTitle(fig), hstr, 

672 self.__saveTempFigure(fig))]) 

673 img_data.extend([('', '', self.__saveTempFigure(x)) 

674 for x in figs[1:]]) 

675 

676 if 'velocity_spectra' in tdict: 

677 fig = self.__createMaxAmpFigure(tid, trcname) 

678 img_data.append((self.__getFigureTitle(fig), '', 

679 self.__saveTempFigure(fig))) 

680 

681 fig = self.__createSpectraFigure(tid, 

682 'velocity_spectra') 

683 img_data.append((self.__getFigureTitle(fig), '', 

684 self.__saveTempFigure(fig))) 

685 if self.output_format == 'pdf': 

686 src_str = self.__formatLatexString( 

687 self.sources[src_id].dump()) 

688 sen_str = self.__formatLatexString( 

689 self.sensors[sen_id].dump()) 

690 else: 

691 src_str = self.sources[src_id].dump() 

692 sen_str = self.sensors[sen_id].dump() 

693 chapters.append([ttl, src_str, sen_str, img_data]) 

694 if self.output_format == 'pdf': 

695 self.__createOutputDoc( 

696 [[self.__formatLatexString(self.store_id), 

697 self.__formatLatexString(self.phase_ratio_string), 

698 artefacts]], chapters) 

699 elif self.output_format == 'html': 

700 self.__createOutputDoc( 

701 [[self.store_id, self.phase_ratio_string, artefacts]], 

702 chapters) 

703 

704 @classmethod 

705 def createComparisonOutputDoc(cls, gft1, gft2, 

706 *trc_ids, **kwargs): 

707 # only valid kwargs is 'together' 

708 if 'together' in kwargs: 

709 together = kwargs['together'] 

710 else: 

711 together = True 

712 

713 if len(trc_ids) == 0: 

714 trc_ids = gft1.traces.keys() 

715 

716 tname = '{0}-{1}{2}{3}'.format( 

717 gft1.store_id, gft2.store_id, '_together_' if together else '_', 

718 gft1.getFrequencyString(None)) 

719 gft1.pdf_name = tname 

720 gft2.pdf_name = tname 

721 

722 trcs1 = gft1.traces 

723 trcs2 = gft2.traces 

724 

725 art1 = gft1.__getArtefactPageInfo(trc_ids) 

726 art2 = gft2.__getArtefactPageInfo(trc_ids, gft1.store_id) 

727 chapters = [] 

728 if gft1.plot_everything: 

729 for tid in trc_ids: 

730 if tid not in trcs1 or tid not in trcs2: 

731 continue 

732 src_id, sen_id = tid.split('|') 

733 ttl = '{0}, {1}'.format( 

734 gft1.getSourceString(src_id), gft1.getSensorString(sen_id)) 

735 tdict1 = trcs1[tid] 

736 tdict2 = trcs2[tid] 

737 img_data = [] 

738 href = r'\hypertarget{${trc_id}|${str_id}|${type}}{}' 

739 trcname = 'displacement_traces' 

740 tstr = 'displacement_spectra' 

741 if trcname in tdict1 and trcname in tdict2: 

742 hstr = Template(href).substitute( 

743 trc_id=tid, str_id=gft1.store_id, type='d') 

744 figs = cls.__createComparisonTraceFigures( 

745 gft1, gft2, tid, trcname, together) 

746 fig = figs[0] 

747 img_data.extend([(gft1.__getFigureTitle(fig), hstr, 

748 gft1.__saveTempFigure(fig))]) 

749 img_data.extend([('', '', gft1.__saveTempFigure(x)) 

750 for x in figs[1:]]) 

751 

752 if tstr in tdict1 and tstr in tdict2: 

753 fig = gft1.__createComparisonMaxAmpFigure( 

754 gft1, gft2, tid, trcname, together) 

755 img_data.append((gft1.__getFigureTitle(fig), '', 

756 gft1.__saveTempFigure(fig))) 

757 

758 trcname = tstr 

759 if not together: 

760 if trcname in tdict1 and trcname in tdict2: 

761 fig = cls.__createComparissonSpectraFigure( 

762 gft1, gft2, tid, trcname) 

763 img_data.append((gft1.__getFigureTitle(fig), '', 

764 gft1.__saveTempFigure(fig))) 

765 

766 trcname = 'velocity_traces' 

767 tstr = 'velocity_spectra' 

768 if gft1.plot_velocity and gft2.plot_velocity and \ 

769 trcname in tdict1 and trcname in tdict2: 

770 

771 hstr = Template(href).substitute( 

772 trc_id=tid, str_id=gft1.store_id, type='v') 

773 figs = cls.__createComparisonTraceFigures( 

774 gft1, gft2, tid, trcname, together) 

775 fig = figs[0] 

776 img_data.extend([(gft1.__getFigureTitle(fig), hstr, 

777 gft1.__saveTempFigure(fig))]) 

778 img_data.extend([('', '', gft1.__saveTempFigure(x)) 

779 for x in figs[1:]]) 

780 

781 if tstr in tdict1 and tstr in tdict2: 

782 fig = gft1.__createComparisonMaxAmpFigure( 

783 gft1, gft2, tid, trcname, together) 

784 img_data.append((gft1.__getFigureTitle(fig), '', 

785 gft1.__saveTempFigure(fig))) 

786 

787 if not together: 

788 if tstr in tdict1 and tstr in tdict2: 

789 fig = cls.__createComparissonSpectraFigure( 

790 gft1, gft2, tid, tstr) 

791 img_data.append((gft1.__getFigureTitle(fig), '', 

792 gft1.__saveTempFigure(fig))) 

793 

794 if gft1.output_format == 'pdf': 

795 src_str = gft1.__formatLatexString( 

796 gft1.sources[src_id].dump()) 

797 sen_str = gft1.__formatLatexString( 

798 gft1.sensors[sen_id].dump()) 

799 else: 

800 src_str = gft1.sources[src_id].dump() 

801 sen_str = gft1.sensors[sen_id].dump() 

802 

803 chapters.append([ttl, src_str, sen_str, img_data]) 

804 

805 if gft1.output_format == 'pdf': 

806 gft1.__createOutputDoc( 

807 [[gft1.__formatLatexString(gft1.store_id), 

808 gft1.__formatLatexString(gft1.phase_ratio_string), art1], 

809 [gft2.__formatLatexString(gft2.store_id), 

810 gft2.__formatLatexString(gft2.phase_ratio_string), art2]], 

811 chapters, gft2=gft2, together=together) 

812 elif gft1.output_format == 'html': 

813 gft1.__createOutputDoc( 

814 [[gft1.store_id, gft1.phase_ratio_string, art1], 

815 [gft2.store_id, gft2.phase_ratio_string, art2]], 

816 chapters, gft2=gft2, together=together) 

817 

818 def __getArtefactPageInfo(self, trc_ids, str_id=None): 

819 is_tex = self.output_format == 'pdf' 

820 if is_tex: 

821 sp = r'\ '*6 

822 else: 

823 sp = r' '*6 

824 

825 ratio_tol = 0.4 

826 artefacts = [] 

827 # list of [<trace string>, <ratio text color>, <ratio string>] 

828 if str_id is None: 

829 str_id = self.store_id 

830 for tid in trc_ids: 

831 if tid not in self.traces: 

832 continue 

833 src_id, sen_id = tid.split('|') 

834 tdict = self.traces[tid] 

835 ttl_str = r'%s, %s' % (self.getSourceString(src_id), 

836 self.getSensorString(sen_id)) 

837 if self.output_format == 'pdf': 

838 ttl_str = r'\textbr{%s}' % ttl_str 

839 

840 artefacts.append([ttl_str, 'black', '']) 

841 

842 chCode = self.sensors[sen_id].codes[3] 

843 ttl_str = r'{0}{1} traces ({2})'.format( 

844 sp, 'Displacement', chCode) 

845 ratio = tdict['displacement_ratio'] 

846 color = ('red' if ratio == 0. or ratio > ratio_tol else 'black') 

847 rat_str = '{0:0.3f}'.format(ratio) 

848 if is_tex: 

849 tex = r'\hyperlink{${tid}|${str_id}|${type}}{${title}}' 

850 ttl_str = Template(tex).substitute( 

851 tid=tid, str_id=str_id, type='d', title=ttl_str) 

852 

853 artefacts.append([ttl_str, color, rat_str]) 

854 

855 if self.plot_velocity and 'velocity_traces' in tdict: 

856 ttl_str = r'{0}{1} traces ({2})'.format(sp, 'Velocity', chCode) 

857 ratio = tdict['velocity_ratio'] 

858 color = ('red' if ratio == 0. or ratio > ratio_tol 

859 else 'black') 

860 rat_str = '{0:0.3f}'.format(ratio) 

861 if is_tex: 

862 tex = r'\hyperlink{${tid}|${str_id}|${type}}{${title}}' 

863 ttl_str = Template(tex).substitute( 

864 tid=tid, str_id=str_id, type='v', title=ttl_str) 

865 

866 artefacts.append([ttl_str, color, rat_str]) 

867 artefacts.append(['', 'black', '']) 

868 return artefacts 

869 

870 def __createTraceFigures(self, trc_id, trc_type): 

871 figs = [] 

872 for i in self.__scalelist: 

873 fig = self.__setupFigure(trc_type, 1) 

874 ax, ax2 = fig.axes 

875 zerotrc = self.__drawTraceData(trc_id, trc_type, ax, ax2, i, 

876 (0.01, 0.01)) 

877 self.__drawBeachball(trc_id, ax) 

878 figs.append(fig) 

879 if zerotrc: 

880 break 

881 return figs 

882 

883 @staticmethod 

884 def __createComparisonTraceFigures( 

885 gfTest1, gfTest2, trc_id, trc_type, together): 

886 

887 tdict = gfTest1.traces[trc_id] 

888 sc1 = tdict['{0}_scale'.format(trc_type.split('_')[0])] 

889 

890 tdict = gfTest2.traces[trc_id] 

891 sc2 = tdict['{0}_scale'.format(trc_type.split('_')[0])] 

892 

893 absmax = (sc1 + sc2) / 2. 

894 figs = [] 

895 for i in gfTest1.__scalelist: 

896 if together: 

897 fig = gfTest1.__setupFigure(trc_type, 1) 

898 ax, ax2 = fig.axes 

899 zerotrc1 = gfTest1.__drawTraceData( 

900 trc_id, trc_type, ax, ax2, i, (0.01, 0.01), 

901 showphases=False, absmax=absmax) 

902 zerotrc2 = gfTest2.__drawTraceData( 

903 trc_id, trc_type, ax, ax2, i, (0.92, 0.01), 

904 color=plot.mpl_color('scarletred2'), hor_ali='right', 

905 showphases=False, absmax=absmax) 

906 gfTest1.__drawBeachball(trc_id, ax) 

907 else: 

908 fig = gfTest1.__setupFigure(trc_type, 2) 

909 ax, ax2, ax3, ax4 = fig.axes 

910 zerotrc1 = gfTest1.__drawTraceData( 

911 trc_id, trc_type, ax, ax2, i, (0.01, 0.01), absmax=absmax) 

912 gfTest1.__drawBeachball(trc_id, ax) 

913 zerotrc2 = gfTest2.__drawTraceData( 

914 trc_id, trc_type, ax3, ax4, i, (0.98, 0.01), 

915 hor_ali='right', absmax=absmax) 

916 gfTest2.__drawBeachball(trc_id, ax3) 

917 

918 figs.append(fig) 

919 if zerotrc1 and zerotrc2: 

920 break 

921 return figs 

922 

923 def __drawTraceData(self, trc_id, trc_type, lfax, rtax, yfac, anno_pt, 

924 color='black', hor_ali='left', showphases=True, 

925 absmax=None): 

926 new_axis = lfax.get_ylim() == (0.0, 1.0) 

927 tdict = self.traces[trc_id] 

928 trcs = tdict[trc_type] 

929 phdict = None 

930 if showphases and 'phase_arrivals' in tdict: 

931 phdict = tdict['phase_arrivals'] 

932 if len(phdict) == 0: 

933 phdict = None 

934 

935 times = trace.minmaxtime(trcs, key=lambda trc: None)[None] 

936 times = np.linspace(times[0], times[1], 10) 

937 diff = (times[1] - times[0]) / 10. 

938 

939 src_id, sen_id = trc_id.split('|') 

940 chCode = self.sensors[sen_id].codes[3] 

941 dists = self.__getSensorDistances(src_id, sen_id) 

942 

943 ysc = dists[1] - dists[0] 

944 

945 zerotrc = False 

946 if absmax is None: 

947 absmax = tdict['{0}_scale'.format(trc_type.split('_')[0])] 

948 if absmax == 0.: 

949 absmax = 1. 

950 zerotrc = True 

951 

952 ysc2 = ysc / absmax * yfac 

953 

954 for i in range(len(self.sensors[sen_id].sensors)): 

955 trc = trcs[i] 

956 dist = dists[i] 

957 lfax.plot(trc.get_xdata(), 

958 (dist + trc.get_ydata() * ysc2) * cake.m2d, '-', 

959 color=color) 

960 

961 if phdict is not None: 

962 for i, pid in enumerate(phdict): 

963 ts, ds = phdict[pid] 

964 ds = [d * cake.m2d for d in ds] 

965 lfax.plot(ts, ds, 

966 label='{0}-wave'.format(pid), marker='o', 

967 markersize=3, 

968 color=plot.to01(plot.graph_colors[i % 7])) 

969 

970 lfax.legend(loc='lower right', shadow=False, prop={'size': 10.}) 

971 

972 xmin = times[0] - diff 

973 xmax = times[-1] + diff 

974 dmin = dists[0] 

975 dmax = dists[-1] 

976 ymin = (dmin - ysc * 3.) 

977 ymax = (dmax + ysc * 3.) 

978 if new_axis: 

979 lfax.set_xlim(xmin, xmax) 

980 lfax.set_ylim(ymin * cake.m2d, ymax * cake.m2d) 

981 rtax.set_ylim(ymin / 1e3, ymax / 1e3) 

982 else: 

983 xlims = lfax.get_xlim() 

984 xlims = (min(xmin, xlims[0]), max(xmax, xlims[1])) 

985 lfax.set_xlim(xlims) 

986 

987 lfax.set_title('{0} traces ({1}), {2}'.format( 

988 trc_type.split('_')[0].title(), chCode, 

989 self.getFrequencyString(trc_id)), y=1.015) 

990 lfax.annotate(self.getAnnotateString(src_id, sen_id, ysc2), 

991 xy=anno_pt, fontsize=(self.__notesize), 

992 xycoords='figure fraction', ha=hor_ali, color=color) 

993 if zerotrc: 

994 lfax.annotate('Zero amplitudes!\nSpectra will not\nbe plotted.', 

995 xy=(0.001, 0.), fontsize=25, alpha=0.75, 

996 rotation=0., xycoords='axes fraction', 

997 color=plot.mpl_color('aluminium4')) 

998 

999 self.__drawGrid(lfax) 

1000 return zerotrc 

1001 

1002 def __createSpectraFigure(self, trc_id, spc_id): 

1003 fig = self.__setupFigure('spectra', 1) 

1004 ax = fig.axes[0] 

1005 self.__drawSpectraData(trc_id, ax, (0.01, 0.01), spc_id) 

1006 return fig 

1007 

1008 @staticmethod 

1009 def __createComparissonSpectraFigure(gfTest1, gfTest2, trc_id, spc_id): 

1010 fig = gfTest1.__setupFigure('spectra', 2) 

1011 ax, ax2 = fig.axes 

1012 gfTest1.__drawSpectraData(trc_id, ax, (0.01, 0.01), spc_id, 

1013 show_cbar=False) 

1014 gfTest2.__drawSpectraData(trc_id, ax2, (0.98, 0.01), spc_id, 

1015 hor_ali='right') 

1016 # evenly space the axes 

1017 ax1, ax2, ax3 = fig.axes 

1018 pos1 = ax1.get_position() 

1019 pos2 = ax2.get_position() 

1020 fac = (pos2.x0 - pos1.x1) / 2. 

1021 mid = (pos2.x1 + pos1.x0) / 2. 

1022 wid = mid - fac - pos1.x0 

1023 ax1.set_position([pos1.x0, pos1.y0, wid, pos1.height]) 

1024 ax2.set_position([mid + fac, pos2.y0, wid, pos2.height]) 

1025 return fig 

1026 

1027 def __drawSpectraData(self, trc_id, ax, anno_pt, spc_id, 

1028 hor_ali='left', show_cbar=True): 

1029 tdict = self.traces[trc_id] 

1030 spcs = tdict[spc_id] 

1031 clrs = np.linspace(0., 1., len(spcs)) 

1032 cmap = cm.jet 

1033 src_id, sen_id = trc_id.split('|') 

1034 chCode = self.sensors[sen_id].codes[3] 

1035 dists = self.__getSensorDistances(src_id, sen_id) 

1036 for idx, data in enumerate(spcs): 

1037 f, v = data 

1038 ax.plot(f, abs(v), color=cmap(clrs[idx]), 

1039 marker='x') 

1040 

1041 ax.set_xscale('log') 

1042 ax.set_yscale('log') 

1043 ax.set_ylabel('Amplitude [{0}]'.format( 

1044 'm' if spc_id.startswith('displacement') else 'm/s')) 

1045 ax.set_xlabel('Frequency [Hz]') 

1046 ax.set_title('{0} spectra for ({1})'. 

1047 format(spc_id.split('_')[0].title(), chCode), y=1.015) 

1048 ax.annotate(self.getAnnotateString(src_id, sen_id), 

1049 xy=anno_pt, ha=hor_ali, fontsize=(self.__notesize), 

1050 xycoords='figure fraction') 

1051 if show_cbar: 

1052 tmin = min(dists) / 1e3 

1053 tmax = max(dists) / 1e3 

1054 sm = plt.cm.ScalarMappable( 

1055 cmap=cmap, norm=plt.Normalize(vmin=tmin, vmax=tmax)) 

1056 sm.set_array(np.linspace(tmin, tmax, 

1057 self.sensors[sen_id].sensor_count)) 

1058 cbar = plt.colorbar(sm, shrink=0.95) 

1059 cbar.ax.set_ylabel('Sensor distance [km]') 

1060 

1061 def __createMaxAmpFigure(self, trc_id, trc_typ): 

1062 fig = self.__setupFigure('maxamp_{0}'.format(trc_typ), 1) 

1063 ax, ax2 = fig.axes 

1064 self.__drawMaxAmpData(trc_id, ax, ax2, (0.01, 0.01), trc_typ) 

1065 return fig 

1066 

1067 @staticmethod 

1068 def __createComparisonMaxAmpFigure(gfTest1, gfTest2, trc_id, trc_typ, 

1069 together): 

1070 if together: 

1071 fig = gfTest1.__setupFigure('maxamp_{0}'.format(trc_typ), 1) 

1072 ax, ax2 = fig.axes 

1073 gfTest1.__drawMaxAmpData(trc_id, ax, ax2, (0.01, 0.01), trc_typ) 

1074 gfTest2.__drawMaxAmpData(trc_id, ax, ax2, (0.92, 0.01), trc_typ, 

1075 color=plot.mpl_color('scarletred2'), 

1076 hor_ali='right') 

1077 else: 

1078 fig = gfTest1.__setupFigure('maxamp_{0}'.format(trc_typ), 2) 

1079 ax, ax2, ax3, ax4 = fig.axes 

1080 gfTest1.__drawMaxAmpData(trc_id, ax, ax2, (0.01, 0.01), trc_typ) 

1081 gfTest2.__drawMaxAmpData(trc_id, ax3, ax4, (0.98, 0.01), trc_typ, 

1082 hor_ali='right') 

1083 return fig 

1084 

1085 def __drawMaxAmpData(self, trc_id, btmax, topax, anno_pt, trc_typ, 

1086 color='black', hor_ali='left'): 

1087 new_axis = btmax.get_ylim() == (0.0, 1.0) 

1088 src_id, sen_id = trc_id.split('|') 

1089 chCode = self.sensors[sen_id].codes[3] 

1090 trcs = self.traces[trc_id][trc_typ] 

1091 dists = [x / 1e3 for x in self.__getSensorDistances(src_id, sen_id)] 

1092 before, after = self.__tracesMinMaxs(trcs, src_id, sen_id) 

1093 amps = [max(abs(x[0]), abs(x[1])) for x in after] 

1094 

1095 btmax.plot(dists, amps, '-x', color=color) 

1096 xlims = [x * cake.m2d * 1e3 for x in btmax.get_xlim()] 

1097 ylims = (min(amps) * 1e-1, max(amps) * 1e1) 

1098 if new_axis: 

1099 topax.set_xlim(xlims) 

1100 btmax.set_ylim(ylims) 

1101 else: 

1102 tlims = btmax.get_ylim() 

1103 ylims = (min(ylims[0], tlims[0]), max(ylims[1], tlims[1])) 

1104 btmax.set_ylim(ylims) 

1105 topax.set_ylim(ylims) 

1106 btmax.set_title('{0} amplitudes for ({1})'. 

1107 format(trc_typ.split('_')[0].title(), chCode), y=1.08) 

1108 btmax.annotate(self.getAnnotateString(src_id, sen_id), 

1109 xy=anno_pt, ha=hor_ali, color=color, 

1110 fontsize=(self.__notesize), xycoords='figure fraction') 

1111 self.__drawGrid(btmax) 

1112 

1113 def __getModelProperty(self, prop): 

1114 mod = self.store.config.earthmodel_1d 

1115 depths = mod.profile('z') / 1e3 

1116 if '/' in prop: 

1117 pros = prop.split('/') 

1118 data = mod.profile(pros[0]) / mod.profile(pros[1]) 

1119 else: 

1120 data = mod.profile(prop) 

1121 

1122 if prop in ['vp', 'vs', 'rho']: 

1123 data /= 1e3 

1124 

1125 return data, depths 

1126 

1127 def __createVelocityProfile(self, gft2=None): 

1128 opts = ['vp', 'vs', 'vp/vs', 'rho', 'qp', 'qs', 'qp/qs'] 

1129 fig = self.__setupFigure('profile', 4, rows=2) 

1130 axs = fig.axes 

1131 for i, opt in enumerate(opts): 

1132 ax = axs[i] 

1133 if gft2 is None: 

1134 data, depths = self.__getModelProperty(opt) 

1135 ax.plot(data, depths) 

1136 else: 

1137 data, depths = self.__getModelProperty(opt) 

1138 ax.plot(data, depths, linewidth=3, 

1139 linestyle='--', alpha=0.8, label=self.store_id, 

1140 color=plot.mpl_color('aluminium4')) 

1141 data, depths = gft2.__getModelProperty(opt) 

1142 ax.plot(data, depths, linewidth=1, 

1143 color=plot.mpl_color('scarletred2'), 

1144 label=gft2.store_id) 

1145 

1146 if opt in ['vp', 'vs']: 

1147 ex = ' [km/s]' 

1148 elif opt == 'rho': 

1149 ex = ' [kg/m^3]' 

1150 else: 

1151 ex = '' 

1152 ax.set_xlabel(opt + ex) 

1153 ax.xaxis.set_label_coords(0.5, -0.13) 

1154 ax.invert_yaxis() 

1155 pos = ax.get_position() 

1156 if i == 0 or i == 4: 

1157 ax.set_ylabel('Depth [km]') 

1158 if i > 1: 

1159 for j in ax.xaxis.get_ticklabels()[1::2]: 

1160 j.set_visible(False) 

1161 if (i // 4) == 0: 

1162 y = pos.y0 * 1.025 

1163 ax.set_position([pos.x0, y, pos.width, 

1164 pos.height - (y - pos.y0)]) 

1165 else: 

1166 y = pos.height * 0.975 

1167 ax.set_position([pos.x0, pos.y0, pos.width, y]) 

1168 

1169 if gft2 is None: 

1170 ttl = 'Earth model plots: {0}'.format(self.store_id) 

1171 else: 

1172 ttl = 'Earth model plots: {0} & {1}'.format(self.store_id, 

1173 gft2.store_id) 

1174 ax.legend(bbox_to_anchor=(1.25, 1.), loc=2) 

1175 

1176 ax.annotate(ttl, xy=(0.5, 0.95), fontsize=(self.__notesize * 2), 

1177 xycoords='figure fraction', ha='center', va='top') 

1178 return self.__saveTempFigure(fig) 

1179 

1180 @staticmethod 

1181 def __setupFigure(fig_type, cols, rows=1): 

1182 fontsize = 10. 

1183 figsize = plot.mpl_papersize('a4', 'landscape') 

1184 plot.mpl_init(fontsize=fontsize) 

1185 

1186 fig = plt.figure(figsize=figsize) 

1187 labelpos = plot.mpl_margins(fig, w=7., h=6., units=fontsize, 

1188 wspace=7., nw=cols, nh=rows) 

1189 for cnt in range(1, (cols * rows) + 1): 

1190 if fig_type == 'profile' and cnt >= 8: 

1191 continue 

1192 ax = fig.add_subplot(rows, cols, cnt) 

1193 labelpos(ax, 2., 2.25) 

1194 if fig_type.startswith('maxamp'): 

1195 if fig_type.split('_')[1] == 'displacement': 

1196 ax.set_ylabel('Max. Amplitude [m]') 

1197 else: 

1198 ax.set_ylabel('Max. Amplitude [m/s]') 

1199 ax.set_yscale('log') 

1200 ax2 = ax.twiny() 

1201 ax.set_xlabel('Distance [km]') 

1202 ax2.set_xlabel('Distance [deg]') 

1203 elif '_traces' in fig_type: 

1204 ax.set_xlabel('Time [s]') 

1205 ax2 = ax.twinx() 

1206 ax2.set_ylabel('Distance [km]') 

1207 ax.set_ylabel('Distance [deg]') 

1208 return fig 

1209 

1210 @staticmethod 

1211 def __drawGrid(ax, major=True, minor=True): 

1212 if major: 

1213 ax.grid( 

1214 visible=True, 

1215 which='major', 

1216 c='grey', 

1217 linestyle='-', 

1218 alpha=.45) 

1219 

1220 if minor: 

1221 ax.minorticks_on() 

1222 ax.grid( 

1223 visible=True, 

1224 which='minor', 

1225 c='grey', 

1226 linestyle=':', 

1227 alpha=.8) 

1228 

1229 @staticmethod 

1230 def __getFigureTitle(fig): 

1231 for ax in fig.axes: 

1232 ttl = ax.get_title() 

1233 if ttl == '': 

1234 continue 

1235 return ttl 

1236 return '' 

1237 

1238 def __drawBeachball(self, trc_id, ax): 

1239 src_id, sen_id = trc_id.split('|') 

1240 src = self.sources[src_id] 

1241 mt = src.pyrocko_moment_tensor() 

1242 

1243 sz = 20. 

1244 szpt = sz / 72. 

1245 bbx = ax.get_xlim()[0] 

1246 bby = ax.get_ylim()[1] 

1247 move_trans = transforms.ScaledTranslation( 

1248 szpt, -szpt, ax.figure.dpi_scale_trans) 

1249 inv_trans = ax.transData.inverted() 

1250 bbx, bby = inv_trans.transform(move_trans.transform( 

1251 ax.transData.transform((bbx, bby)))) 

1252 beachball.plot_beachball_mpl(mt, ax, beachball_type='full', 

1253 size=sz, position=(bbx, bby)) 

1254 

1255 def __saveTempFigure(self, fig): 

1256 if self.output_format == 'pdf': 

1257 fname = NamedTemporaryFile( 

1258 prefix='gft_', suffix='.pdf', dir=self.temp_dir, delete=False) 

1259 fname.close() 

1260 fig.savefig(fname.name, transparent=True) 

1261 out = fname.name 

1262 elif self.output_format == 'html': 

1263 imgdata = io.BytesIO() 

1264 fig.savefig(imgdata, format='png') 

1265 imgdata.seek(0) 

1266 out = base64.b64encode(imgdata.read()) 

1267 plt.close(fig) 

1268 return out 

1269 

1270 def __tracesMinMaxs(self, trcs, src_id, sen_id): 

1271 # return the min/max amplitudes before and after the arrival of the 

1272 # self.phase_ratio_string phase found in the traces, 

1273 # plus the maximum ratio between the max absolute value 

1274 before = [] 

1275 after = [] 

1276 

1277 tbrk = None 

1278 src = self.sources[src_id] 

1279 sensors = self.sensors[sen_id].sensors 

1280 for i, trc in enumerate(trcs): 

1281 if self.phases == '': 

1282 tbrk = None 

1283 else: 

1284 tbrk = self.store.t(self.phase_ratio_string, src, sensors[i]) 

1285 times = trc.get_xdata() 

1286 data = trc.get_ydata() 

1287 tmina = None 

1288 tmaxa = None 

1289 tminb = None 

1290 tmaxb = None 

1291 for idx, time in enumerate(times): 

1292 val = data[idx] 

1293 if time < tbrk or tbrk is None: 

1294 if tminb is None or tminb > val: 

1295 tminb = val 

1296 if tmaxb is None or tmaxb < val: 

1297 tmaxb = val 

1298 else: 

1299 if tmina is None or tmina > val: 

1300 tmina = val 

1301 if tmaxa is None or tmaxa < val: 

1302 tmaxa = val 

1303 if tminb is None: 

1304 tminb = 0. 

1305 if tmaxb is None: 

1306 tmaxb = 0. 

1307 before.append((tminb, tmaxb)) 

1308 if tbrk is None: 

1309 after.append((tminb, tmaxb)) 

1310 else: 

1311 after.append((tmina, tmaxa)) 

1312 

1313 return before, after 

1314 

1315 def __tracesMinMax(self, trcs, src_id, sen_id): 

1316 # return the min/max amplitudes before and after the arrival of the 

1317 # self.phase_ratio_string phase found in the traces, 

1318 # plus the maximum ratio between the max absolute value 

1319 before, after = self.__tracesMinMaxs(trcs, src_id, sen_id) 

1320 mina = min([x[0] for x in after]) 

1321 maxa = max([x[1] for x in after]) 

1322 minb = min([x[0] for x in before]) 

1323 maxb = max([x[1] for x in before]) 

1324 ratios = list(map(lambda b, a: 0. if a[0] == a[1] == 0. else 

1325 max(abs(b[0]), abs(b[1]))/max(abs(a[0]), abs(a[1])), 

1326 before, after)) 

1327 

1328 return mina, maxa, minb, maxb, max(ratios) 

1329 

1330 def __getSensorDistances(self, src_id, sen_id): 

1331 src = self.sources[src_id] 

1332 return [src.distance_to(sen) for sen in self.sensors[sen_id].sensors] 

1333 

1334 @classmethod 

1335 def __createStandardSetups(cls, store_dir, **kwargs): 

1336 tdict = cls.__get_valid_arguments(kwargs) 

1337 

1338 gft = cls(store_dir, **tdict) 

1339 

1340 if 'source_depth' in kwargs: 

1341 depth = kwargs['source_depth'] 

1342 else: 

1343 depth = None 

1344 src1 = gft.createSource('DC', depth, 0., 90., 0., **kwargs) 

1345 # src2 = gft.createSource('DC', depth, -90., 90., -90., **kwargs) 

1346 src3 = gft.createSource('DC', depth, 45., 90., 180., **kwargs) 

1347 # src4 = gft.createSource('Explosion', depth, 0., 90., **kwargs) 

1348 

1349 SensorArray.validate_args(kwargs) 

1350 sen1 = gft.createSensors(strike=0., codes=('', 'STA', '', 'R'), 

1351 azimuth=0., dip=0., **kwargs) 

1352 sen2 = gft.createSensors(strike=0., codes=('', 'STA', '', 'T'), 

1353 azimuth=90., dip=0., **kwargs) 

1354 sen3 = gft.createSensors(strike=0., codes=('', 'STA', '', 'Z'), 

1355 azimuth=0., dip=-90., **kwargs) 

1356 

1357 gft.trace_configs = [[src3, sen1], [src1, sen2], [src3, sen3]] 

1358 # gft.trace_configs = [[src3, sen1]] 

1359 # gft.trace_configs = [[src3, sen1], [src1, sen2], [src3, sen3], 

1360 # [src4, 'all']] 

1361 # gft.trace_configs = [[src4, 'all']] 

1362 # gft.trace_configs = [['all', 'all']] 

1363 gft.__applyTypicalProcedures() 

1364 return gft 

1365 

1366 def __applyTypicalProcedures(self): 

1367 if self.trace_configs is None: 

1368 self.createDisplacementTraces() 

1369 else: 

1370 for src, sen in self.trace_configs: 

1371 self.createDisplacementTraces(src, sen) 

1372 if self.plot_velocity: 

1373 if self.trace_configs is None: 

1374 self.createVelocityTraces() 

1375 else: 

1376 for src, sen in self.trace_configs: 

1377 self.createVelocityTraces('{0}|{1}'.format(src, sen)) 

1378 self.applyFrequencyFilters() 

1379 self.getPhaseArrivals() 

1380 

1381 def __update(self, **kwargs): 

1382 for k in kwargs: 

1383 temp = kwargs[k] 

1384 if temp is not None: 

1385 setattr(self, k, temp) 

1386 

1387 @classmethod 

1388 def runStandardCheck( 

1389 cls, store_dir, source_depth=None, 

1390 lowpass_frequency=None, highpass_frequency=None, 

1391 rel_lowpass_frequency=(1. / 110), rel_highpass_frequency=(1. / 16), 

1392 distance_min=None, distance_max=None, sensor_count=50, 

1393 filter_order=4, pdf_dir=None, plot_velocity=False, 

1394 plot_everything=True, output_format='pdf'): 

1395 

1396 args = locals() 

1397 del args['cls'] 

1398 gft = cls.__createStandardSetups(**args) 

1399 gft.createOutputDoc() 

1400 gft.__printMessage() 

1401 return gft 

1402 

1403 @classmethod 

1404 def runComparissonStandardCheck( 

1405 cls, store_dir1, store_dir2, distance_min, distance_max, 

1406 together=True, source_depth=None, output_format='pdf', 

1407 lowpass_frequency=None, highpass_frequency=None, 

1408 rel_lowpass_frequency=(1. / 110), rel_highpass_frequency=(1. / 16), 

1409 filter_order=4, pdf_dir=None, plot_velocity=False, 

1410 plot_everything=True, sensor_count=50): 

1411 

1412 args = locals() 

1413 del args['cls'] 

1414 args['change_dists'] = False 

1415 gft1 = cls.__createStandardSetups(store_dir1, **args) 

1416 if 'source_depth' not in args or args['source_depth'] is None: 

1417 args['source_depth'] = gft1.sources[list(gft1.sources.keys())[0]]\ 

1418 .depth 

1419 

1420 gft2 = cls.__createStandardSetups(store_dir2, **args) 

1421 

1422 cls.createComparisonOutputDoc(gft1, gft2, together=together) 

1423 gft1.__printMessage() 

1424 gft2.__printMessage() 

1425 return gft1, gft2 

1426 

1427 @staticmethod 

1428 def createDocumentFromFile(filename, allowed=1, **kwargs): 

1429 with open(filename, 'r') as f: 

1430 fstr = f.read() 

1431 gfts = [] 

1432 st = None 

1433 end = -1 

1434 tstr = '--- !{0}.GreensFunctionTest'.format(guts_prefix) 

1435 try: 

1436 while True: 

1437 st = fstr.index(tstr, end + 1) 

1438 if st is None: 

1439 break 

1440 if st > 0 and fstr[st - 1:st] != '\n': 

1441 end = st 

1442 st = None 

1443 continue 

1444 end = fstr.index('\n{0}'.format(tstr), st) 

1445 gfts.append(fstr[st:end]) 

1446 except ValueError: 

1447 if st is not None: 

1448 gfts.append(fstr[st:]) 

1449 

1450 cnt = len(gfts) 

1451 if cnt == 0: 

1452 raise TypeError('error: GreensFunctionTest: non-valid config' 

1453 ' file passed: {0}'.format(filename)) 

1454 elif cnt > allowed: 

1455 raise TypeError('error: GreensFunctionTest: to many configs in' 

1456 ' file passed: {0}'.format(filename)) 

1457 elif cnt < allowed: 

1458 raise TypeError('error: GreensFunctionTest: not enough configs in' 

1459 ' file passed: {0}'.format(filename)) 

1460 

1461 if cnt == allowed == 1: 

1462 gft = load(stream=gfts[0]) 

1463 gft.__update(**kwargs) 

1464 gft.__applyTypicalProcedures() 

1465 gft.createOutputDoc() 

1466 gft.__printMessage() 

1467 return gft 

1468 

1469 if cnt == allowed == 2: 

1470 gft = load(stream=gfts[0]) 

1471 gft.__update(**kwargs) 

1472 gft.__applyTypicalProcedures() 

1473 

1474 gft2 = load(stream=gfts[1]) 

1475 gft2.__update(**kwargs) 

1476 gft2.__applyTypicalProcedures() 

1477 

1478 gft.createComparisonOutputDoc(gft, gft2, **kwargs) 

1479 gft.__printMessage() 

1480 gft2.__printMessage() 

1481 return gft, gft2