1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
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
18import numpy as np
19from matplotlib import pyplot as plt
20from matplotlib import cm, transforms
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
29from jinja2 import Environment, PackageLoader
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'))
46class FomostoReportError(Exception):
47 pass
50class GreensFunctionError(FomostoReportError):
51 pass
54class FilterFrequencyError(GreensFunctionError):
56 def __init__(self, frequency):
57 Exception.__init__(self)
58 self.frequency = frequency
60 def __str__(self):
61 return 'Cannot set {0} frequency to both an absolute and' \
62 ' revlative value.'.format(self.frequency)
65class SourceError(GreensFunctionError):
67 def __init__(self, typ):
68 Exception.__init__(self)
69 self.type = typ
71 def __str__(self):
72 return 'Source type not currently supported: {0}'.format(self.type)
75class SensorArray(Target):
77 distance_min = Float.T()
78 distance_max = Float.T()
79 strike = Float.T()
80 sensor_count = Int.T(default=50)
82 __target_allowed = dir(Target)
83 __sensorarray_allowed = ['distance_min', 'distance_max', 'strike',
84 'sensor_count', 'change_dists']
86 def __init__(self, **kwargs):
87 Object.__init__(self, **kwargs)
89 self.validate_args(kwargs)
90 self.__build(**kwargs)
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)
100 for i in ks:
101 del kwargs[i]
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)
110 for i in ks:
111 del kwargs[i]
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]
126class GreensFunctionTest(Object):
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
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)
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
162 def __init__(self, store_dir, **kwargs):
163 Object.__init__(self, store_dir=store_dir, **kwargs)
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')
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)
177 if self.rel_lowpass_frequency is not None:
178 self.lowpass_frequency = self.store.config.sample_rate * \
179 self.rel_lowpass_frequency
181 if self.rel_highpass_frequency is not None:
182 self.highpass_frequency = self.store.config.sample_rate * \
183 self.rel_highpass_frequency
185 self.rel_lowpass_frequency = None
186 self.rel_highpass_frequency = None
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)
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)
203 self.traces = OrderedDict()
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()
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
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'
225 def __printMessage(self):
226 if self.message is not None:
227 print(self.message)
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)
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
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
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
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'
277 if tid not in self.traces:
278 return ''
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'
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
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))
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)
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)
353 self.sources[tstr] = src
354 self.src_ids.append(tstr)
355 return tstr
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))
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))
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))
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)
403 self.sensors[tstr] = sen_ar
404 self.sen_ids.append(tstr)
405 return tstr
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))
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))
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)
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)
497 if len(times) > 0:
498 phdict[pid] = (times, dists)
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
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
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
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
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))
546 def __createOutputDoc(self, artefacts, chapters,
547 gft2=None, together=False):
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
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()
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
575 info = [(fstr_id, config, tpath, img_ttl)]
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)
585 config = gft2.store.config.dump()
586 if together:
587 tpath = ''
588 else:
589 tpath = gft2.__createVelocityProfile()
591 info += [(fstr_id, config, tpath, fstr_id)]
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')))
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?')
616 self.cleanup()
617 if gft2 is not None:
618 gft2.cleanup()
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)
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:]])
655 if 'displacement_spectra' in tdict:
656 fig = self.__createMaxAmpFigure(tid, trcname)
657 img_data.append((self.__getFigureTitle(fig), '',
658 self.__saveTempFigure(fig)))
660 fig = self.__createSpectraFigure(tid,
661 'displacement_spectra')
662 img_data.append((self.__getFigureTitle(fig), '',
663 self.__saveTempFigure(fig)))
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:]])
676 if 'velocity_spectra' in tdict:
677 fig = self.__createMaxAmpFigure(tid, trcname)
678 img_data.append((self.__getFigureTitle(fig), '',
679 self.__saveTempFigure(fig)))
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)
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
713 if len(trc_ids) == 0:
714 trc_ids = gft1.traces.keys()
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
722 trcs1 = gft1.traces
723 trcs2 = gft2.traces
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:]])
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)))
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)))
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:
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:]])
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)))
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)))
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()
803 chapters.append([ttl, src_str, sen_str, img_data])
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)
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
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
840 artefacts.append([ttl_str, 'black', ''])
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)
853 artefacts.append([ttl_str, color, rat_str])
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)
866 artefacts.append([ttl_str, color, rat_str])
867 artefacts.append(['', 'black', ''])
868 return artefacts
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
883 @staticmethod
884 def __createComparisonTraceFigures(
885 gfTest1, gfTest2, trc_id, trc_type, together):
887 tdict = gfTest1.traces[trc_id]
888 sc1 = tdict['{0}_scale'.format(trc_type.split('_')[0])]
890 tdict = gfTest2.traces[trc_id]
891 sc2 = tdict['{0}_scale'.format(trc_type.split('_')[0])]
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)
918 figs.append(fig)
919 if zerotrc1 and zerotrc2:
920 break
921 return figs
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
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.
939 src_id, sen_id = trc_id.split('|')
940 chCode = self.sensors[sen_id].codes[3]
941 dists = self.__getSensorDistances(src_id, sen_id)
943 ysc = dists[1] - dists[0]
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
952 ysc2 = ysc / absmax * yfac
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)
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]))
970 lfax.legend(loc='lower right', shadow=False, prop={'size': 10.})
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)
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'))
999 self.__drawGrid(lfax)
1000 return zerotrc
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
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
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')
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]')
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
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
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]
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)
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)
1122 if prop in ['vp', 'vs', 'rho']:
1123 data /= 1e3
1125 return data, depths
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)
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])
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)
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)
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)
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
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)
1220 if minor:
1221 ax.minorticks_on()
1222 ax.grid(
1223 visible=True,
1224 which='minor',
1225 c='grey',
1226 linestyle=':',
1227 alpha=.8)
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 ''
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()
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))
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
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 = []
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))
1313 return before, after
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))
1328 return mina, maxa, minb, maxb, max(ratios)
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]
1334 @classmethod
1335 def __createStandardSetups(cls, store_dir, **kwargs):
1336 tdict = cls.__get_valid_arguments(kwargs)
1338 gft = cls(store_dir, **tdict)
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)
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)
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
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()
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)
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'):
1396 args = locals()
1397 del args['cls']
1398 gft = cls.__createStandardSetups(**args)
1399 gft.createOutputDoc()
1400 gft.__printMessage()
1401 return gft
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):
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
1420 gft2 = cls.__createStandardSetups(store_dir2, **args)
1422 cls.createComparisonOutputDoc(gft1, gft2, together=together)
1423 gft1.__printMessage()
1424 gft2.__printMessage()
1425 return gft1, gft2
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:])
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))
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
1469 if cnt == allowed == 2:
1470 gft = load(stream=gfts[0])
1471 gft.__update(**kwargs)
1472 gft.__applyTypicalProcedures()
1474 gft2 = load(stream=gfts[1])
1475 gft2.__update(**kwargs)
1476 gft2.__applyTypicalProcedures()
1478 gft.createComparisonOutputDoc(gft, gft2, **kwargs)
1479 gft.__printMessage()
1480 gft2.__printMessage()
1481 return gft, gft2