1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, print_function, division
7import os
8import sys
9import shutil
10import subprocess
11import re
12import io
13import base64
14import datetime
15from collections import OrderedDict
16from tempfile import NamedTemporaryFile, mkdtemp
17from string import Template
19import numpy as np
20from matplotlib import pyplot as plt
21from matplotlib import cm, transforms
23from pyrocko import gf, trace, cake, util, plot
24from pyrocko.response import DifferentiationResponse
25from pyrocko.plot import beachball
26from pyrocko.guts import load, Object, String, List, Float, Int, Bool, Dict
27from pyrocko.gf import Source, Target
28from pyrocko.gf.meta import OutOfBounds
30from jinja2 import Environment, PackageLoader
32guts_prefix = 'gft'
33ex_path = os.path.dirname(os.path.abspath(sys.argv[0]))
34ja_latex_env = Environment(block_start_string='\\BLOCK{',
35 block_end_string='}',
36 variable_start_string='\\VAR{',
37 variable_end_string='}',
38 comment_start_string='\\%{',
39 comment_end_string='}',
40 line_statement_prefix='\\-',
41 line_comment_prefix='%%',
42 trim_blocks=True,
43 loader=PackageLoader('pyrocko',
44 'data/fomosto_report'))
47class FomostoReportError(Exception):
48 pass
51class GreensFunctionError(FomostoReportError):
52 pass
55class FilterFrequencyError(GreensFunctionError):
57 def __init__(self, frequency):
58 Exception.__init__(self)
59 self.frequency = frequency
61 def __str__(self):
62 return 'Cannot set {0} frequency to both an absolute and' \
63 ' revlative value.'.format(self.frequency)
66class SourceError(GreensFunctionError):
68 def __init__(self, typ):
69 Exception.__init__(self)
70 self.type = typ
72 def __str__(self):
73 return 'Source type not currently supported: {0}'.format(self.type)
76class SensorArray(Target):
78 distance_min = Float.T()
79 distance_max = Float.T()
80 strike = Float.T()
81 sensor_count = Int.T(default=50)
83 __target_allowed = dir(Target)
84 __sensorarray_allowed = ['distance_min', 'distance_max', 'strike',
85 'sensor_count', 'change_dists']
87 def __init__(self, **kwargs):
88 Object.__init__(self, **kwargs)
90 self.validate_args(kwargs)
91 self.__build(**kwargs)
93 @classmethod
94 def validate_args(cls, kwargs):
95 ks = []
96 for k in kwargs:
97 if not (k in cls.__target_allowed or
98 k in cls.__sensorarray_allowed):
99 ks.append(k)
101 for i in ks:
102 del kwargs[i]
104 @classmethod
105 def __validate_args_target(cls, kwargs):
106 ks = []
107 for k in kwargs:
108 if k not in cls.__target_allowed:
109 ks.append(k)
111 for i in ks:
112 del kwargs[i]
114 def __build(self, **kwargs):
115 if 'quantity' not in kwargs:
116 kwargs['quantity'] = 'displacement'
117 dists = np.linspace(self.distance_min, self.distance_max,
118 self.sensor_count)
119 self.__validate_args_target(kwargs)
120 self.sensors = [
121 Target(north_shift=float(np.cos(np.radians(self.strike)) * dist),
122 east_shift=float(np.sin(np.radians(self.strike)) * dist),
123 **kwargs)
124 for dist in dists]
127class GreensFunctionTest(Object):
129 __valid_properties = ['store_dir', 'pdf_dir', 'output_format',
130 'lowpass_frequency', 'rel_lowpass_frequency',
131 'highpass_frequency', 'rel_highpass_frequency',
132 'filter_order',
133 'plot_velocity', 'plot_everything']
134 __notesize = 7.45
135 __scalelist = [1, 5, 9.5, 19, 29]
136 __has_phase = True
138 store_dir = String.T()
139 pdf_dir = String.T(default=os.path.expanduser('~'))
140 output_format = String.T(optional=True, default='pdf')
141 lowpass_frequency = Float.T(optional=True)
142 highpass_frequency = Float.T(optional=True)
143 rel_lowpass_frequency = Float.T(optional=True)
144 rel_highpass_frequency = Float.T(optional=True)
145 filter_order = Int.T(default=4)
146 sensor_count = Int.T(default=50)
147 plot_velocity = Bool.T(default=False)
148 plot_everything = Bool.T(default=True)
149 src_ids = List.T(String.T())
150 sen_ids = List.T(String.T())
151 sources = Dict.T(String.T(), Source.T())
152 sensors = Dict.T(String.T(), SensorArray.T())
153 trace_configs = List.T(List.T(String.T()), optional=True)
155 @classmethod
156 def __get_valid_arguments(cls, args):
157 tdict = {}
158 for k in args:
159 if k in cls.__valid_properties:
160 tdict[k] = args[k]
161 return tdict
163 def __init__(self, store_dir, **kwargs):
164 Object.__init__(self, store_dir=store_dir, **kwargs)
166 if self.lowpass_frequency and self.rel_lowpass_frequency:
167 raise FilterFrequencyError('lowpass')
168 if self.highpass_frequency and self.rel_highpass_frequency:
169 raise FilterFrequencyError('highpass')
171 if self.store_dir[-1] != '/':
172 self.store_dir += '/'
173 self.engine = gf.LocalEngine(store_dirs=[self.store_dir])
174 ids = self.engine.get_store_ids()
175 self.store_id = ids[0]
176 self.store = self.engine.get_store(self.store_id)
178 if self.rel_lowpass_frequency is not None:
179 self.lowpass_frequency = self.store.config.sample_rate * \
180 self.rel_lowpass_frequency
182 if self.rel_highpass_frequency is not None:
183 self.highpass_frequency = self.store.config.sample_rate * \
184 self.rel_highpass_frequency
186 self.rel_lowpass_frequency = None
187 self.rel_highpass_frequency = None
189 self.phases = '|'.join([x.id
190 for x in self.store.config.tabulated_phases])
191 if self.phases == '':
192 self.phase_ratio_string = r'\color{red}' \
193 '{warning: store has no tabulated phases}'
194 elif 'begin' in self.phases:
195 self.phase_ratio_string = 'begin'
196 else:
197 self.phase_ratio_string = 'first({0})'.format(self.phases)
199 if len(self.src_ids) == 0 and len(self.sources) > 0:
200 self.src_ids = sorted(self.sources)
201 if len(self.sen_ids) == 0 and len(self.sensors) > 0:
202 self.sen_ids = sorted(self.sensors)
204 self.traces = OrderedDict()
206 if self.pdf_dir is None:
207 self.pdf_dir = os.path.expanduser('~')
208 if self.pdf_dir[-1] != '/':
209 self.pdf_dir += '/'
210 self.pdf_name = '{0}_{1}'.format(
211 self.store_id, self.getFrequencyString(None))
212 util.setup_logging()
214 self.temp_dir = mkdtemp(prefix='gft_')
215 self.message = None
216 self.changed_depth = False
217 self.changed_dist_min = False
218 self.changed_dist_max = False
220 def __addToMessage(self, msg):
221 if self.message is None:
222 self.message = '\nmessage(s) for store {0}:\n'. \
223 format(self.store_id)
224 self.message += msg + '\n'
226 def __printMessage(self):
227 if self.message is not None:
228 print(self.message)
230 def cleanup(self):
231 shutil.rmtree(self.temp_dir)
232 for i in ['.aux', '.log', '.out', '.toc']:
233 path = '{0}{1}{2}'.format(self.pdf_dir, self.pdf_name, i)
234 if os.path.exists(path):
235 os.remove(path)
237 def getSourceString(self, sid):
238 if sid not in self.sources:
239 return ''
240 src = self.sources[sid]
241 if isinstance(src, gf.DCSource):
242 typ = 'Double Couple'
243 sstr = 'Source Type: {0}, Strike: {1:g}, Dip: {2:g}, Rake: {3:g}' \
244 .format(typ, src.strike, src.dip, src.rake)
245 elif isinstance(src, gf.RectangularExplosionSource):
246 typ = 'Explosion'
247 sstr = 'Source Type: {0}, Strike: {1:g}, Dip: {2:g}'.format(
248 typ, src.strike, src.dip)
249 else:
250 typ = '{0}'.format(type(src)).split('.')[-1].split("'")[0]
251 sstr = 'Source Type: {0}, see config page'.format(typ)
252 return sstr
254 def getSensorString(self, sid):
255 if sid not in self.sensors:
256 return ''
257 sen_ar = self.sensors[sid]
258 sensors = sen_ar.sensors
260 targ = sensors[0]
261 senstr = 'Sensor Strike: {0:g}, Azimuth: {1:g}, Dip: {2:g}'. \
262 format(sen_ar.strike, targ.azimuth, targ.dip)
263 return senstr
265 def getFrequencyString(self, tid):
266 lpf = self.lowpass_frequency
267 hpf = self.highpass_frequency
268 if tid is None:
269 if lpf and hpf:
270 return '{0:.4g}-{1:.4g}Hz'.format(lpf, hpf)
271 elif lpf:
272 return 'lowpass {0:.4g}Hz'.format(lpf)
273 elif hpf:
274 return 'highpass {0:.4g}Hz'.format(hpf)
275 else:
276 return 'unfiltered'
278 if tid not in self.traces:
279 return ''
281 tdict = self.traces[tid]
282 lpfa = tdict['lowpass_applied']
283 hpfa = tdict['highpass_applied']
284 if lpf and hpf and lpfa and hpfa:
285 return '{0:.4g}-{1:.4g} [Hz]'.format(lpf, hpf)
286 elif lpf and lpfa:
287 return 'lowpass frequency: {0:.4g} [Hz]'.format(lpf)
288 elif hpf and hpfa:
289 return 'highpass frequency: {0:.4g} [Hz]'.format(hpf)
290 else:
291 return 'unfiltered'
293 def getAnnotateString(self, src_id, sen_id, fac=None):
294 tstr = 'Green\'s Function: {0}\n{1}\n{2}'.format(
295 self.store_id, self.getSourceString(src_id),
296 self.getSensorString(sen_id))
297 if fac is not None:
298 tstr += '\nAmplitude Scale Factor: {0:0.4g} [microns/km]'. \
299 format(fac * 1e-9 / 1e3)
300 return tstr
302 def createSource(self, typ, depth, strike, dip, rake=None, lat=0., lon=0.,
303 north_shift=0., east_shift=0., change_dists=True,
304 **kwargs):
305 stCfg = self.store.config
306 stMin = stCfg.source_depth_min
307 stMax = stCfg.source_depth_max
308 if depth is None:
309 depth = (stMax + stMin) / 2.
310 elif depth < stMin or depth > stMax:
311 if not change_dists:
312 raise OutOfBounds('Source depth.')
313 diff = abs(stMin - depth)
314 if diff < stMax - depth:
315 depth = stMin
316 else:
317 depth = stMax
318 if not self.changed_depth:
319 self.changed_depth = True
320 self.__addToMessage(
321 'Source depth out of bounds. Changed to: {0:g}m.'.
322 format(depth))
324 if typ.upper() in ('EX' 'EXPLOSION'):
325 src = gf.RectangularExplosionSource(
326 lat=lat,
327 lon=lon,
328 north_shift=north_shift,
329 east_shift=east_shift,
330 depth=depth,
331 strike=strike,
332 dip=dip)
333 elif typ.upper() in ('DC' 'DOUBLE COUPLE'):
334 src = gf.DCSource(
335 lat=lat,
336 lon=lon,
337 north_shift=north_shift,
338 east_shift=east_shift,
339 depth=depth,
340 strike=strike,
341 dip=dip,
342 rake=rake)
343 else:
344 raise SourceError(typ)
346 srcstr = '{0}.{1:.2g}.{2:.2g}.{3:.2g}'.format(
347 typ[0:2], depth, strike, dip)
348 tstr = srcstr
349 tint = 96
350 while tstr in self.sources:
351 tint += 1
352 tstr = srcstr + chr(tint)
354 self.sources[tstr] = src
355 self.src_ids.append(tstr)
356 return tstr
358 def createSensors(self, distance_min=None, distance_max=None,
359 change_dists=True, **kwargs):
360 stCfg = self.store.config
361 stMin = stCfg.distance_min
362 stMax = stCfg.distance_max
363 if distance_min is None:
364 distance_min = stMin
365 elif distance_min < stMin or distance_min > stMax:
366 if not change_dists:
367 raise OutOfBounds('Minimum sensor distance.')
368 distance_min = stMin
369 if not self.changed_dist_min:
370 self.changed_dist_min = True
371 self.__addToMessage(
372 'Sensor minimum distance out of bounds. Changed to allowed'
373 ' minimum {0:g}m.'.format(stMin))
375 if distance_max is None:
376 distance_max = stMax
377 elif (distance_max > stMax or distance_max < stMin):
378 if not change_dists:
379 raise OutOfBounds('Maximum sensor distance')
380 distance_max = stMax
381 if not self.changed_dist_max:
382 self.changed_dist_max = True
383 self.__addToMessage(
384 'Sensor maximum distance out of bounds. Changed to allowed'
385 ' maximum {0:g}m.'.format(stMax))
387 if change_dists and distance_min == distance_max:
388 distance_min = stMin
389 distance_max = stMax
390 self.__addToMessage(
391 'Sensor minimum and maximum distance equal. Changed to'
392 ' allowed minimum {0:g}m and maximum {1:g}m.'.format(stMin,
393 stMax))
395 sen_ar = SensorArray(distance_min=distance_min,
396 distance_max=distance_max, **kwargs)
397 senstr = '{0}.{1:.2g}'.format('.'.join(sen_ar.codes), sen_ar.strike)
398 tstr = senstr
399 tint = 96
400 while tstr in self.sensors:
401 tint += 1
402 tstr = senstr + chr(tint)
404 self.sensors[tstr] = sen_ar
405 self.sen_ids.append(tstr)
406 return tstr
408 def createDisplacementTraces(self, src_id='all', sen_id='all'):
409 for sr_id in self.src_ids:
410 if sr_id not in self.sources:
411 continue
412 if not(src_id == 'all' or sr_id == src_id):
413 continue
414 for sn_id in self.sen_ids:
415 if sn_id not in self.sensors:
416 continue
417 if not(sen_id == 'all' or sn_id == sen_id):
418 continue
419 try:
420 response = self.engine.process(
421 self.sources[sr_id],
422 self.sensors[sn_id].sensors)
423 trcs = response.pyrocko_traces()
424 tstr = '{0}|{1}'.format(sr_id, sn_id)
425 if tstr not in self.traces:
426 self.traces[tstr] = {}
427 tdict = self.traces[tstr]
428 tdict['displacement_traces'] = trcs
429 mina, maxa, minb, maxb, ratio = \
430 self.__tracesMinMax(trcs, sr_id, sn_id)
431 if ratio != 0.:
432 tdict['displacement_spectra'] = [
433 trc.spectrum() for trc in trcs]
434 tdict['lowpass_applied'] = False
435 tdict['highpass_applied'] = False
436 tdict['displacement_ratio'] = ratio
437 tdict['displacement_scale'] = max(abs(mina), abs(maxa))
438 except IndexError:
439 self.__addToMessage(
440 'warning: IndexError: no traces created for'
441 ' source-sensor combination: {0} - {1}. Try increasing'
442 ' the sensor minimum distance.'.format(
443 sr_id, sn_id))
445 def createVelocityTraces(self, trc_id='all'):
446 for tid in self.traces:
447 ids = trc_id.split('|')
448 if len(ids) == 2:
449 src_id, sen_id = ids
450 else:
451 src_id = sen_id = None
452 if not(trc_id == 'all' or
453 ((src_id == 'all' or src_id in self.sources) and
454 (sen_id == 'all' or sen_id in self.sensors))):
455 continue
456 tdict = self.traces[tid]
457 if 'displacement_traces' not in tdict:
458 continue
459 trcs = [trc.transfer(
460 transfer_function=DifferentiationResponse())
461 for trc in tdict['displacement_traces']]
462 tdict['velocity_traces'] = trcs
463 src_id, sen_id = tid.split('|')
464 mina, maxa, minb, maxb, ratio = \
465 self.__tracesMinMax(trcs, src_id, sen_id)
466 if ratio != 0.:
467 tdict['velocity_spectra'] = [
468 trc.spectrum() for trc in trcs]
469 tdict['velocity_ratio'] = ratio
470 tdict['velocity_scale'] = max(abs(mina), abs(maxa))
472 def getPhaseArrivals(self, trc_id='all', phase_ids='all'):
473 for tid in self.traces:
474 if trc_id == 'all' or trc_id == tid:
475 self.__setPhaseArrivals(tid, phase_ids)
477 def __setPhaseArrivals(self, trc_id, phase_ids='all'):
478 st = self.store
479 src_id, sen_id = trc_id.split('|')
480 src = self.sources[src_id]
481 tdict = self.traces[trc_id]
482 if 'phase_arrivals' not in tdict:
483 tdict['phase_arrivals'] = {}
484 phdict = tdict['phase_arrivals']
485 for pid in self.phases.split('|'):
486 if pid == '':
487 continue
488 if phase_ids == 'all' or pid == phase_ids:
489 dists = []
490 times = []
491 for sen in self.sensors[sen_id].sensors:
492 time = st.t(pid, src, sen)
493 if time is None:
494 continue
495 dists.append(src.distance_to(sen))
496 times.append(time)
498 if len(times) > 0:
499 phdict[pid] = (times, dists)
501 def applyFrequencyFilters(self, trc_id='all'):
502 dispname = 'displacement_traces'
503 velname = 'velocity_traces'
504 for tid in self.traces:
505 if trc_id == 'all' or trc_id == tid:
506 tdict = self.traces[tid]
507 if self.lowpass_frequency:
508 if dispname in tdict:
509 for trc in tdict[dispname]:
510 trc.lowpass(self.filter_order,
511 self.lowpass_frequency, demean=False)
512 tdict['lowpass_applied'] = True
514 if velname in tdict:
515 for trc in tdict[velname]:
516 trc.lowpass(self.filter_order,
517 self.lowpass_frequency, demean=False)
518 tdict['lowpass_applied'] = True
520 if self.highpass_frequency:
521 if dispname in tdict:
522 for trc in tdict[dispname]:
523 trc.highpass(self.filter_order,
524 self.highpass_frequency, demean=False)
525 tdict['highpass_applied'] = True
527 if velname in tdict:
528 for trc in tdict[velname]:
529 trc.highpass(self.filter_order,
530 self.highpass_frequency, demean=False)
531 tdict['highpass_applied'] = True
533 sr_id, sn_id = tid.split('|')
534 if dispname in tdict:
535 trcs = tdict[dispname]
536 mina, maxa, minb, maxb, ratio = \
537 self.__tracesMinMax(trcs, sr_id, sn_id)
538 tdict['displacement_ratio'] = ratio
539 tdict['displacement_scale'] = max(abs(mina), abs(maxa))
540 if velname in tdict:
541 trcs = tdict[velname]
542 mina, maxa, minb, maxb, ratio = \
543 self.__tracesMinMax(trcs, sr_id, sn_id)
544 tdict['velocity_ratio'] = ratio
545 tdict['velocity_scale'] = max(abs(mina), abs(maxa))
547 def __createOutputDoc(self, artefacts, chapters,
548 gft2=None, together=False):
550 str_id = self.store_id
551 is_tex = self.output_format == 'pdf'
552 if is_tex:
553 file_type = 'tex'
554 dir = self.temp_dir
555 fstr_id = self.__formatLatexString(str_id)
556 else:
557 file_type = self.output_format
558 dir = self.pdf_dir
559 fstr_id = str_id
561 temp = ja_latex_env.get_template('gfreport.{0}'.format(file_type))
562 out = '{0}/{1}.{2}'.format(dir, self.pdf_name, file_type)
563 config = self.store.config.dump()
565 if together:
566 tpath = self.__createVelocityProfile(gft2)
567 if is_tex:
568 img_ttl = r'{0} \& {1}'.format(
569 fstr_id, gft2.__formatLatexString(gft2.store_id))
570 else:
571 img_ttl = r'{0} & {1}'.format(str_id, gft2.store_id)
572 else:
573 tpath = self.__createVelocityProfile()
574 img_ttl = fstr_id
576 info = [(fstr_id, config, tpath, img_ttl)]
578 if gft2 is not None:
579 if is_tex:
580 fstr_id = self.__formatLatexString(gft2.store_id)
581 str_id = r'{0} \& {1}'.format(str_id, gft2.store_id)
582 else:
583 fstr_id = gft2.store_id
584 str_id = r'{0} & {1}'.format(str_id, gft2.store_id)
586 config = gft2.store.config.dump()
587 if together:
588 tpath = ''
589 else:
590 tpath = gft2.__createVelocityProfile()
592 info += [(fstr_id, config, tpath, fstr_id)]
594 with open(out, 'w') as f:
595 if is_tex:
596 f.write(temp.render(
597 rpt_id=self.__formatLatexString(str_id), str_info=info,
598 artefacts=artefacts, chapters=chapters,
599 config=config, headings='headings'))
600 elif file_type == 'html':
601 f.write(temp.render(
602 rpt_id=str_id, str_info=info, artefacts=artefacts,
603 chapters=chapters, config=config,
604 date=datetime.datetime.now().strftime('%B %d, %Y')))
606 if is_tex:
607 pro_call = ['pdflatex', '-output-directory', self.pdf_dir, out,
608 '-interaction', 'nonstop']
609 try:
610 subprocess.call(pro_call)
611 subprocess.call(pro_call)
612 subprocess.call(pro_call)
613 except OSError:
614 raise FomostoReportError(
615 'Cannot run "pdflatex" executable. Is it installed?')
617 self.cleanup()
618 if gft2 is not None:
619 gft2.cleanup()
621 @staticmethod
622 def __formatLatexString(string):
623 rep = {'_': r'\_', r'\n': r'\\', '|': r'\textbar '}
624 rep = {re.escape(k): v for k, v in rep.items()}
625 pat = re.compile('|'.join(rep.keys()))
626 return pat.sub(lambda m: rep[re.escape(m.group(0))], string)
628 def createOutputDoc(self, *trc_ids):
629 trcs = self.traces
630 if len(trc_ids) == 0:
631 trc_ids = trcs.keys()
632 artefacts = self.__getArtefactPageInfo(trc_ids)
633 chapters = []
634 if self.plot_everything:
635 for tid in trc_ids:
636 if tid not in trcs:
637 continue
638 src_id, sen_id = tid.split('|')
639 ttl = '{0}, {1}'.format(
640 self.getSourceString(src_id),
641 self.getSensorString(sen_id))
642 tdict = trcs[tid]
643 img_data = []
644 href = r'\hypertarget{${trc_id}|${str_id}|${type}}{}'
645 trcname = 'displacement_traces'
646 if trcname in tdict:
647 hstr = Template(href).substitute(
648 trc_id=tid, str_id=self.store_id, type='d')
649 figs = self.__createTraceFigures(tid, trcname)
650 fig = figs[0]
651 img_data.extend([(self.__getFigureTitle(fig), hstr,
652 self.__saveTempFigure(fig))])
653 img_data.extend([('', '', self.__saveTempFigure(x))
654 for x in figs[1:]])
656 if 'displacement_spectra' in tdict:
657 fig = self.__createMaxAmpFigure(tid, trcname)
658 img_data.append((self.__getFigureTitle(fig), '',
659 self.__saveTempFigure(fig)))
661 fig = self.__createSpectraFigure(tid,
662 'displacement_spectra')
663 img_data.append((self.__getFigureTitle(fig), '',
664 self.__saveTempFigure(fig)))
666 trcname = 'velocity_traces'
667 if self.plot_velocity and trcname in tdict:
668 hstr = Template(href).substitute(
669 trc_id=tid, str_id=self.store_id, type='v')
670 figs = self.__createTraceFigures(tid, trcname)
671 fig = figs[0]
672 img_data.extend([(self.__getFigureTitle(fig), hstr,
673 self.__saveTempFigure(fig))])
674 img_data.extend([('', '', self.__saveTempFigure(x))
675 for x in figs[1:]])
677 if 'velocity_spectra' in tdict:
678 fig = self.__createMaxAmpFigure(tid, trcname)
679 img_data.append((self.__getFigureTitle(fig), '',
680 self.__saveTempFigure(fig)))
682 fig = self.__createSpectraFigure(tid,
683 'velocity_spectra')
684 img_data.append((self.__getFigureTitle(fig), '',
685 self.__saveTempFigure(fig)))
686 if self.output_format == 'pdf':
687 src_str = self.__formatLatexString(
688 self.sources[src_id].dump())
689 sen_str = self.__formatLatexString(
690 self.sensors[sen_id].dump())
691 else:
692 src_str = self.sources[src_id].dump()
693 sen_str = self.sensors[sen_id].dump()
694 chapters.append([ttl, src_str, sen_str, img_data])
695 if self.output_format == 'pdf':
696 self.__createOutputDoc(
697 [[self.__formatLatexString(self.store_id),
698 self.__formatLatexString(self.phase_ratio_string),
699 artefacts]], chapters)
700 elif self.output_format == 'html':
701 self.__createOutputDoc(
702 [[self.store_id, self.phase_ratio_string, artefacts]],
703 chapters)
705 @classmethod
706 def createComparisonOutputDoc(cls, gft1, gft2,
707 *trc_ids, **kwargs):
708 # only valid kwargs is 'together'
709 if 'together' in kwargs:
710 together = kwargs['together']
711 else:
712 together = True
714 if len(trc_ids) == 0:
715 trc_ids = gft1.traces.keys()
717 tname = '{0}-{1}{2}{3}'.format(
718 gft1.store_id, gft2.store_id, '_together_' if together else '_',
719 gft1.getFrequencyString(None))
720 gft1.pdf_name = tname
721 gft2.pdf_name = tname
723 trcs1 = gft1.traces
724 trcs2 = gft2.traces
726 art1 = gft1.__getArtefactPageInfo(trc_ids)
727 art2 = gft2.__getArtefactPageInfo(trc_ids, gft1.store_id)
728 chapters = []
729 if gft1.plot_everything:
730 for tid in trc_ids:
731 if tid not in trcs1 or tid not in trcs2:
732 continue
733 src_id, sen_id = tid.split('|')
734 ttl = '{0}, {1}'.format(
735 gft1.getSourceString(src_id), gft1.getSensorString(sen_id))
736 tdict1 = trcs1[tid]
737 tdict2 = trcs2[tid]
738 img_data = []
739 href = r'\hypertarget{${trc_id}|${str_id}|${type}}{}'
740 trcname = 'displacement_traces'
741 tstr = 'displacement_spectra'
742 if trcname in tdict1 and trcname in tdict2:
743 hstr = Template(href).substitute(
744 trc_id=tid, str_id=gft1.store_id, type='d')
745 figs = cls.__createComparisonTraceFigures(
746 gft1, gft2, tid, trcname, together)
747 fig = figs[0]
748 img_data.extend([(gft1.__getFigureTitle(fig), hstr,
749 gft1.__saveTempFigure(fig))])
750 img_data.extend([('', '', gft1.__saveTempFigure(x))
751 for x in figs[1:]])
753 if tstr in tdict1 and tstr in tdict2:
754 fig = gft1.__createComparisonMaxAmpFigure(
755 gft1, gft2, tid, trcname, together)
756 img_data.append((gft1.__getFigureTitle(fig), '',
757 gft1.__saveTempFigure(fig)))
759 trcname = tstr
760 if not together:
761 if trcname in tdict1 and trcname in tdict2:
762 fig = cls.__createComparissonSpectraFigure(
763 gft1, gft2, tid, trcname)
764 img_data.append((gft1.__getFigureTitle(fig), '',
765 gft1.__saveTempFigure(fig)))
767 trcname = 'velocity_traces'
768 tstr = 'velocity_spectra'
769 if gft1.plot_velocity and gft2.plot_velocity and \
770 trcname in tdict1 and trcname in tdict2:
772 hstr = Template(href).substitute(
773 trc_id=tid, str_id=gft1.store_id, type='v')
774 figs = cls.__createComparisonTraceFigures(
775 gft1, gft2, tid, trcname, together)
776 fig = figs[0]
777 img_data.extend([(gft1.__getFigureTitle(fig), hstr,
778 gft1.__saveTempFigure(fig))])
779 img_data.extend([('', '', gft1.__saveTempFigure(x))
780 for x in figs[1:]])
782 if tstr in tdict1 and tstr in tdict2:
783 fig = gft1.__createComparisonMaxAmpFigure(
784 gft1, gft2, tid, trcname, together)
785 img_data.append((gft1.__getFigureTitle(fig), '',
786 gft1.__saveTempFigure(fig)))
788 if not together:
789 if tstr in tdict1 and tstr in tdict2:
790 fig = cls.__createComparissonSpectraFigure(
791 gft1, gft2, tid, tstr)
792 img_data.append((gft1.__getFigureTitle(fig), '',
793 gft1.__saveTempFigure(fig)))
795 if gft1.output_format == 'pdf':
796 src_str = gft1.__formatLatexString(
797 gft1.sources[src_id].dump())
798 sen_str = gft1.__formatLatexString(
799 gft1.sensors[sen_id].dump())
800 else:
801 src_str = gft1.sources[src_id].dump()
802 sen_str = gft1.sensors[sen_id].dump()
804 chapters.append([ttl, src_str, sen_str, img_data])
806 if gft1.output_format == 'pdf':
807 gft1.__createOutputDoc(
808 [[gft1.__formatLatexString(gft1.store_id),
809 gft1.__formatLatexString(gft1.phase_ratio_string), art1],
810 [gft2.__formatLatexString(gft2.store_id),
811 gft2.__formatLatexString(gft2.phase_ratio_string), art2]],
812 chapters, gft2=gft2, together=together)
813 elif gft1.output_format == 'html':
814 gft1.__createOutputDoc(
815 [[gft1.store_id, gft1.phase_ratio_string, art1],
816 [gft2.store_id, gft2.phase_ratio_string, art2]],
817 chapters, gft2=gft2, together=together)
819 def __getArtefactPageInfo(self, trc_ids, str_id=None):
820 is_tex = self.output_format == 'pdf'
821 if is_tex:
822 sp = r'\ '*6
823 else:
824 sp = r' '*6
826 ratio_tol = 0.4
827 artefacts = []
828 # list of [<trace string>, <ratio text color>, <ratio string>]
829 if str_id is None:
830 str_id = self.store_id
831 for tid in trc_ids:
832 if tid not in self.traces:
833 continue
834 src_id, sen_id = tid.split('|')
835 tdict = self.traces[tid]
836 ttl_str = r'%s, %s' % (self.getSourceString(src_id),
837 self.getSensorString(sen_id))
838 if self.output_format == 'pdf':
839 ttl_str = r'\textbr{%s}' % ttl_str
841 artefacts.append([ttl_str, 'black', ''])
843 chCode = self.sensors[sen_id].codes[3]
844 ttl_str = r'{0}{1} traces ({2})'.format(
845 sp, 'Displacement', chCode)
846 ratio = tdict['displacement_ratio']
847 color = ('red' if ratio == 0. or ratio > ratio_tol else 'black')
848 rat_str = '{0:0.3f}'.format(ratio)
849 if is_tex:
850 tex = r'\hyperlink{${tid}|${str_id}|${type}}{${title}}'
851 ttl_str = Template(tex).substitute(
852 tid=tid, str_id=str_id, type='d', title=ttl_str)
854 artefacts.append([ttl_str, color, rat_str])
856 if self.plot_velocity and 'velocity_traces' in tdict:
857 ttl_str = r'{0}{1} traces ({2})'.format(sp, 'Velocity', chCode)
858 ratio = tdict['velocity_ratio']
859 color = ('red' if ratio == 0. or ratio > ratio_tol
860 else 'black')
861 rat_str = '{0:0.3f}'.format(ratio)
862 if is_tex:
863 tex = r'\hyperlink{${tid}|${str_id}|${type}}{${title}}'
864 ttl_str = Template(tex).substitute(
865 tid=tid, str_id=str_id, type='v', title=ttl_str)
867 artefacts.append([ttl_str, color, rat_str])
868 artefacts.append(['', 'black', ''])
869 return artefacts
871 def __createTraceFigures(self, trc_id, trc_type):
872 figs = []
873 for i in self.__scalelist:
874 fig = self.__setupFigure(trc_type, 1)
875 ax, ax2 = fig.axes
876 zerotrc = self.__drawTraceData(trc_id, trc_type, ax, ax2, i,
877 (0.01, 0.01))
878 self.__drawBeachball(trc_id, ax)
879 figs.append(fig)
880 if zerotrc:
881 break
882 return figs
884 @staticmethod
885 def __createComparisonTraceFigures(
886 gfTest1, gfTest2, trc_id, trc_type, together):
888 tdict = gfTest1.traces[trc_id]
889 sc1 = tdict['{0}_scale'.format(trc_type.split('_')[0])]
891 tdict = gfTest2.traces[trc_id]
892 sc2 = tdict['{0}_scale'.format(trc_type.split('_')[0])]
894 absmax = (sc1 + sc2) / 2.
895 figs = []
896 for i in gfTest1.__scalelist:
897 if together:
898 fig = gfTest1.__setupFigure(trc_type, 1)
899 ax, ax2 = fig.axes
900 zerotrc1 = gfTest1.__drawTraceData(
901 trc_id, trc_type, ax, ax2, i, (0.01, 0.01),
902 showphases=False, absmax=absmax)
903 zerotrc2 = gfTest2.__drawTraceData(
904 trc_id, trc_type, ax, ax2, i, (0.92, 0.01),
905 color=plot.mpl_color('scarletred2'), hor_ali='right',
906 showphases=False, absmax=absmax)
907 gfTest1.__drawBeachball(trc_id, ax)
908 else:
909 fig = gfTest1.__setupFigure(trc_type, 2)
910 ax, ax2, ax3, ax4 = fig.axes
911 zerotrc1 = gfTest1.__drawTraceData(
912 trc_id, trc_type, ax, ax2, i, (0.01, 0.01), absmax=absmax)
913 gfTest1.__drawBeachball(trc_id, ax)
914 zerotrc2 = gfTest2.__drawTraceData(
915 trc_id, trc_type, ax3, ax4, i, (0.98, 0.01),
916 hor_ali='right', absmax=absmax)
917 gfTest2.__drawBeachball(trc_id, ax3)
919 figs.append(fig)
920 if zerotrc1 and zerotrc2:
921 break
922 return figs
924 def __drawTraceData(self, trc_id, trc_type, lfax, rtax, yfac, anno_pt,
925 color='black', hor_ali='left', showphases=True,
926 absmax=None):
927 new_axis = lfax.get_ylim() == (0.0, 1.0)
928 tdict = self.traces[trc_id]
929 trcs = tdict[trc_type]
930 phdict = None
931 if showphases and 'phase_arrivals' in tdict:
932 phdict = tdict['phase_arrivals']
933 if len(phdict) == 0:
934 phdict = None
936 times = trace.minmaxtime(trcs, key=lambda trc: None)[None]
937 times = np.linspace(times[0], times[1], 10)
938 diff = (times[1] - times[0]) / 10.
940 src_id, sen_id = trc_id.split('|')
941 chCode = self.sensors[sen_id].codes[3]
942 dists = self.__getSensorDistances(src_id, sen_id)
944 ysc = dists[1] - dists[0]
946 zerotrc = False
947 if absmax is None:
948 absmax = tdict['{0}_scale'.format(trc_type.split('_')[0])]
949 if absmax == 0.:
950 absmax = 1.
951 zerotrc = True
953 ysc2 = ysc / absmax * yfac
955 for i in range(len(self.sensors[sen_id].sensors)):
956 trc = trcs[i]
957 dist = dists[i]
958 lfax.plot(trc.get_xdata(),
959 (dist + trc.get_ydata() * ysc2) * cake.m2d, '-',
960 color=color)
962 if phdict is not None:
963 for i, pid in enumerate(phdict):
964 ts, ds = phdict[pid]
965 ds = [d * cake.m2d for d in ds]
966 lfax.plot(ts, ds,
967 label='{0}-wave'.format(pid), marker='o',
968 markersize=3,
969 color=plot.to01(plot.graph_colors[i % 7]))
971 lfax.legend(loc='lower right', shadow=False, prop={'size': 10.})
973 xmin = times[0] - diff
974 xmax = times[-1] + diff
975 dmin = dists[0]
976 dmax = dists[-1]
977 ymin = (dmin - ysc * 3.)
978 ymax = (dmax + ysc * 3.)
979 if new_axis:
980 lfax.set_xlim(xmin, xmax)
981 lfax.set_ylim(ymin * cake.m2d, ymax * cake.m2d)
982 rtax.set_ylim(ymin / 1e3, ymax / 1e3)
983 else:
984 xlims = lfax.get_xlim()
985 xlims = (min(xmin, xlims[0]), max(xmax, xlims[1]))
986 lfax.set_xlim(xlims)
988 lfax.set_title('{0} traces ({1}), {2}'.format(
989 trc_type.split('_')[0].title(), chCode,
990 self.getFrequencyString(trc_id)), y=1.015)
991 lfax.annotate(self.getAnnotateString(src_id, sen_id, ysc2),
992 xy=anno_pt, fontsize=(self.__notesize),
993 xycoords='figure fraction', ha=hor_ali, color=color)
994 if zerotrc:
995 lfax.annotate("Zero amplitudes!\nSpectra will not\nbe plotted.",
996 xy=(0.001, 0.), fontsize=25, alpha=0.75,
997 rotation=0., xycoords='axes fraction',
998 color=plot.mpl_color('aluminium4'))
1000 self.__drawGrid(lfax)
1001 return zerotrc
1003 def __createSpectraFigure(self, trc_id, spc_id):
1004 fig = self.__setupFigure('spectra', 1)
1005 ax = fig.axes[0]
1006 self.__drawSpectraData(trc_id, ax, (0.01, 0.01), spc_id)
1007 return fig
1009 @staticmethod
1010 def __createComparissonSpectraFigure(gfTest1, gfTest2, trc_id, spc_id):
1011 fig = gfTest1.__setupFigure('spectra', 2)
1012 ax, ax2 = fig.axes
1013 gfTest1.__drawSpectraData(trc_id, ax, (0.01, 0.01), spc_id,
1014 show_cbar=False)
1015 gfTest2.__drawSpectraData(trc_id, ax2, (0.98, 0.01), spc_id,
1016 hor_ali='right')
1017 # evenly space the axes
1018 ax1, ax2, ax3 = fig.axes
1019 pos1 = ax1.get_position()
1020 pos2 = ax2.get_position()
1021 fac = (pos2.x0 - pos1.x1) / 2.
1022 mid = (pos2.x1 + pos1.x0) / 2.
1023 wid = mid - fac - pos1.x0
1024 ax1.set_position([pos1.x0, pos1.y0, wid, pos1.height])
1025 ax2.set_position([mid + fac, pos2.y0, wid, pos2.height])
1026 return fig
1028 def __drawSpectraData(self, trc_id, ax, anno_pt, spc_id,
1029 hor_ali='left', show_cbar=True):
1030 tdict = self.traces[trc_id]
1031 spcs = tdict[spc_id]
1032 clrs = np.linspace(0., 1., len(spcs))
1033 cmap = cm.jet
1034 src_id, sen_id = trc_id.split('|')
1035 chCode = self.sensors[sen_id].codes[3]
1036 dists = self.__getSensorDistances(src_id, sen_id)
1037 for idx, data in enumerate(spcs):
1038 f, v = data
1039 ax.plot(f, abs(v), color=cmap(clrs[idx]),
1040 marker='x')
1042 ax.set_xscale('log')
1043 ax.set_yscale('log')
1044 ax.set_ylabel('Amplitude [{0}]'.format(
1045 'm' if spc_id.startswith('displacement') else 'm/s'))
1046 ax.set_xlabel('Frequency [Hz]')
1047 ax.set_title('{0} spectra for ({1})'.
1048 format(spc_id.split('_')[0].title(), chCode), y=1.015)
1049 ax.annotate(self.getAnnotateString(src_id, sen_id),
1050 xy=anno_pt, ha=hor_ali, fontsize=(self.__notesize),
1051 xycoords='figure fraction')
1052 if show_cbar:
1053 tmin = min(dists) / 1e3
1054 tmax = max(dists) / 1e3
1055 sm = plt.cm.ScalarMappable(
1056 cmap=cmap, norm=plt.Normalize(vmin=tmin, vmax=tmax))
1057 sm.set_array(np.linspace(tmin, tmax,
1058 self.sensors[sen_id].sensor_count))
1059 cbar = plt.colorbar(sm, shrink=0.95)
1060 cbar.ax.set_ylabel('Sensor distance [km]')
1062 def __createMaxAmpFigure(self, trc_id, trc_typ):
1063 fig = self.__setupFigure('maxamp_{0}'.format(trc_typ), 1)
1064 ax, ax2 = fig.axes
1065 self.__drawMaxAmpData(trc_id, ax, ax2, (0.01, 0.01), trc_typ)
1066 return fig
1068 @staticmethod
1069 def __createComparisonMaxAmpFigure(gfTest1, gfTest2, trc_id, trc_typ,
1070 together):
1071 if together:
1072 fig = gfTest1.__setupFigure('maxamp_{0}'.format(trc_typ), 1)
1073 ax, ax2 = fig.axes
1074 gfTest1.__drawMaxAmpData(trc_id, ax, ax2, (0.01, 0.01), trc_typ)
1075 gfTest2.__drawMaxAmpData(trc_id, ax, ax2, (0.92, 0.01), trc_typ,
1076 color=plot.mpl_color('scarletred2'),
1077 hor_ali='right')
1078 else:
1079 fig = gfTest1.__setupFigure('maxamp_{0}'.format(trc_typ), 2)
1080 ax, ax2, ax3, ax4 = fig.axes
1081 gfTest1.__drawMaxAmpData(trc_id, ax, ax2, (0.01, 0.01), trc_typ)
1082 gfTest2.__drawMaxAmpData(trc_id, ax3, ax4, (0.98, 0.01), trc_typ,
1083 hor_ali='right')
1084 return fig
1086 def __drawMaxAmpData(self, trc_id, btmax, topax, anno_pt, trc_typ,
1087 color='black', hor_ali='left'):
1088 new_axis = btmax.get_ylim() == (0.0, 1.0)
1089 src_id, sen_id = trc_id.split('|')
1090 chCode = self.sensors[sen_id].codes[3]
1091 trcs = self.traces[trc_id][trc_typ]
1092 dists = [x / 1e3 for x in self.__getSensorDistances(src_id, sen_id)]
1093 before, after = self.__tracesMinMaxs(trcs, src_id, sen_id)
1094 amps = [max(abs(x[0]), abs(x[1])) for x in after]
1096 btmax.plot(dists, amps, '-x', color=color)
1097 xlims = [x * cake.m2d * 1e3 for x in btmax.get_xlim()]
1098 ylims = (min(amps) * 1e-1, max(amps) * 1e1)
1099 if new_axis:
1100 topax.set_xlim(xlims)
1101 btmax.set_ylim(ylims)
1102 else:
1103 tlims = btmax.get_ylim()
1104 ylims = (min(ylims[0], tlims[0]), max(ylims[1], tlims[1]))
1105 btmax.set_ylim(ylims)
1106 topax.set_ylim(ylims)
1107 btmax.set_title('{0} amplitudes for ({1})'.
1108 format(trc_typ.split('_')[0].title(), chCode), y=1.08)
1109 btmax.annotate(self.getAnnotateString(src_id, sen_id),
1110 xy=anno_pt, ha=hor_ali, color=color,
1111 fontsize=(self.__notesize), xycoords='figure fraction')
1112 self.__drawGrid(btmax)
1114 def __getModelProperty(self, prop):
1115 mod = self.store.config.earthmodel_1d
1116 depths = mod.profile('z') / 1e3
1117 if '/' in prop:
1118 pros = prop.split('/')
1119 data = mod.profile(pros[0]) / mod.profile(pros[1])
1120 else:
1121 data = mod.profile(prop)
1123 if prop in ['vp', 'vs', 'rho']:
1124 data /= 1e3
1126 return data, depths
1128 def __createVelocityProfile(self, gft2=None):
1129 opts = ['vp', 'vs', 'vp/vs', 'rho', 'qp', 'qs', 'qp/qs']
1130 fig = self.__setupFigure('profile', 4, rows=2)
1131 axs = fig.axes
1132 for i, opt in enumerate(opts):
1133 ax = axs[i]
1134 if gft2 is None:
1135 data, depths = self.__getModelProperty(opt)
1136 ax.plot(data, depths)
1137 else:
1138 data, depths = self.__getModelProperty(opt)
1139 ax.plot(data, depths, linewidth=3,
1140 linestyle='--', alpha=0.8, label=self.store_id,
1141 color=plot.mpl_color('aluminium4'))
1142 data, depths = gft2.__getModelProperty(opt)
1143 ax.plot(data, depths, linewidth=1,
1144 color=plot.mpl_color('scarletred2'),
1145 label=gft2.store_id)
1147 if opt in ['vp', 'vs']:
1148 ex = ' [km/s]'
1149 elif opt == 'rho':
1150 ex = ' [kg/m^3]'
1151 else:
1152 ex = ''
1153 ax.set_xlabel(opt + ex)
1154 ax.xaxis.set_label_coords(0.5, -0.13)
1155 ax.invert_yaxis()
1156 pos = ax.get_position()
1157 if i == 0 or i == 4:
1158 ax.set_ylabel('Depth [km]')
1159 if i > 1:
1160 for j in ax.xaxis.get_ticklabels()[1::2]:
1161 j.set_visible(False)
1162 if (i // 4) == 0:
1163 y = pos.y0 * 1.025
1164 ax.set_position([pos.x0, y, pos.width,
1165 pos.height - (y - pos.y0)])
1166 else:
1167 y = pos.height * 0.975
1168 ax.set_position([pos.x0, pos.y0, pos.width, y])
1170 if gft2 is None:
1171 ttl = 'Earth model plots: {0}'.format(self.store_id)
1172 else:
1173 ttl = 'Earth model plots: {0} & {1}'.format(self.store_id,
1174 gft2.store_id)
1175 ax.legend(bbox_to_anchor=(1.25, 1.), loc=2)
1177 ax.annotate(ttl, xy=(0.5, 0.95), fontsize=(self.__notesize * 2),
1178 xycoords='figure fraction', ha='center', va='top')
1179 return self.__saveTempFigure(fig)
1181 @staticmethod
1182 def __setupFigure(fig_type, cols, rows=1):
1183 fontsize = 10.
1184 figsize = plot.mpl_papersize('a4', 'landscape')
1185 plot.mpl_init(fontsize=fontsize)
1187 fig = plt.figure(figsize=figsize)
1188 labelpos = plot.mpl_margins(fig, w=7., h=6., units=fontsize,
1189 wspace=7., nw=cols, nh=rows)
1190 for cnt in range(1, (cols * rows) + 1):
1191 if fig_type == 'profile' and cnt >= 8:
1192 continue
1193 ax = fig.add_subplot(rows, cols, cnt)
1194 labelpos(ax, 2., 2.25)
1195 if fig_type.startswith('maxamp'):
1196 if fig_type.split('_')[1] == 'displacement':
1197 ax.set_ylabel('Max. Amplitude [m]')
1198 else:
1199 ax.set_ylabel('Max. Amplitude [m/s]')
1200 ax.set_yscale('log')
1201 ax2 = ax.twiny()
1202 ax.set_xlabel('Distance [km]')
1203 ax2.set_xlabel('Distance [deg]')
1204 elif '_traces' in fig_type:
1205 ax.set_xlabel('Time [s]')
1206 ax2 = ax.twinx()
1207 ax2.set_ylabel('Distance [km]')
1208 ax.set_ylabel('Distance [deg]')
1209 return fig
1211 @staticmethod
1212 def __drawGrid(ax, major=True, minor=True):
1213 if major:
1214 ax.grid(b=True, which='major', c='grey', linestyle='-', alpha=.45)
1215 if minor:
1216 ax.minorticks_on()
1217 ax.grid(b=True, which='minor', c='grey', linestyle=':', alpha=.8)
1219 @staticmethod
1220 def __getFigureTitle(fig):
1221 for ax in fig.axes:
1222 ttl = ax.get_title()
1223 if ttl == '':
1224 continue
1225 return ttl
1226 return ''
1228 def __drawBeachball(self, trc_id, ax):
1229 src_id, sen_id = trc_id.split('|')
1230 src = self.sources[src_id]
1231 mt = src.pyrocko_moment_tensor()
1233 sz = 20.
1234 szpt = sz / 72.
1235 bbx = ax.get_xlim()[0]
1236 bby = ax.get_ylim()[1]
1237 move_trans = transforms.ScaledTranslation(
1238 szpt, -szpt, ax.figure.dpi_scale_trans)
1239 inv_trans = ax.transData.inverted()
1240 bbx, bby = inv_trans.transform(move_trans.transform(
1241 ax.transData.transform((bbx, bby))))
1242 beachball.plot_beachball_mpl(mt, ax, beachball_type='full',
1243 size=sz, position=(bbx, bby))
1245 def __saveTempFigure(self, fig):
1246 if self.output_format == 'pdf':
1247 fname = NamedTemporaryFile(
1248 prefix='gft_', suffix='.pdf', dir=self.temp_dir, delete=False)
1249 fname.close()
1250 fig.savefig(fname.name, transparent=True)
1251 out = fname.name
1252 elif self.output_format == 'html':
1253 imgdata = io.BytesIO()
1254 fig.savefig(imgdata, format='png')
1255 imgdata.seek(0)
1256 out = base64.b64encode(imgdata.read())
1257 plt.close(fig)
1258 return out
1260 def __tracesMinMaxs(self, trcs, src_id, sen_id):
1261 # return the min/max amplitudes before and after the arrival of the
1262 # self.phase_ratio_string phase found in the traces,
1263 # plus the maximum ratio between the max absolute value
1264 before = []
1265 after = []
1267 tbrk = None
1268 src = self.sources[src_id]
1269 sensors = self.sensors[sen_id].sensors
1270 for i, trc in enumerate(trcs):
1271 if self.phases == '':
1272 tbrk = None
1273 else:
1274 tbrk = self.store.t(self.phase_ratio_string, src, sensors[i])
1275 times = trc.get_xdata()
1276 data = trc.get_ydata()
1277 tmina = None
1278 tmaxa = None
1279 tminb = None
1280 tmaxb = None
1281 for idx, time in enumerate(times):
1282 val = data[idx]
1283 if time < tbrk or tbrk is None:
1284 if tminb is None or tminb > val:
1285 tminb = val
1286 if tmaxb is None or tmaxb < val:
1287 tmaxb = val
1288 else:
1289 if tmina is None or tmina > val:
1290 tmina = val
1291 if tmaxa is None or tmaxa < val:
1292 tmaxa = val
1293 if tminb is None:
1294 tminb = 0.
1295 if tmaxb is None:
1296 tmaxb = 0.
1297 before.append((tminb, tmaxb))
1298 if tbrk is None:
1299 after.append((tminb, tmaxb))
1300 else:
1301 after.append((tmina, tmaxa))
1303 return before, after
1305 def __tracesMinMax(self, trcs, src_id, sen_id):
1306 # return the min/max amplitudes before and after the arrival of the
1307 # self.phase_ratio_string phase found in the traces,
1308 # plus the maximum ratio between the max absolute value
1309 before, after = self.__tracesMinMaxs(trcs, src_id, sen_id)
1310 mina = min([x[0] for x in after])
1311 maxa = max([x[1] for x in after])
1312 minb = min([x[0] for x in before])
1313 maxb = max([x[1] for x in before])
1314 ratios = list(map(lambda b, a: 0. if a[0] == a[1] == 0. else
1315 max(abs(b[0]), abs(b[1]))/max(abs(a[0]), abs(a[1])),
1316 before, after))
1318 return mina, maxa, minb, maxb, max(ratios)
1320 def __getSensorDistances(self, src_id, sen_id):
1321 src = self.sources[src_id]
1322 return [src.distance_to(sen) for sen in self.sensors[sen_id].sensors]
1324 @classmethod
1325 def __createStandardSetups(cls, store_dir, **kwargs):
1326 tdict = cls.__get_valid_arguments(kwargs)
1328 gft = cls(store_dir, **tdict)
1330 if 'source_depth' in kwargs:
1331 depth = kwargs['source_depth']
1332 else:
1333 depth = None
1334 src1 = gft.createSource('DC', depth, 0., 90., 0., **kwargs)
1335 # src2 = gft.createSource('DC', depth, -90., 90., -90., **kwargs)
1336 src3 = gft.createSource('DC', depth, 45., 90., 180., **kwargs)
1337 # src4 = gft.createSource('Explosion', depth, 0., 90., **kwargs)
1339 SensorArray.validate_args(kwargs)
1340 sen1 = gft.createSensors(strike=0., codes=('', 'STA', '', 'R'),
1341 azimuth=0., dip=0., **kwargs)
1342 sen2 = gft.createSensors(strike=0., codes=('', 'STA', '', 'T'),
1343 azimuth=90., dip=0., **kwargs)
1344 sen3 = gft.createSensors(strike=0., codes=('', 'STA', '', 'Z'),
1345 azimuth=0., dip=-90., **kwargs)
1347 gft.trace_configs = [[src3, sen1], [src1, sen2], [src3, sen3]]
1348 # gft.trace_configs = [[src3, sen1]]
1349 # gft.trace_configs = [[src3, sen1], [src1, sen2], [src3, sen3],
1350 # [src4, 'all']]
1351 # gft.trace_configs = [[src4, 'all']]
1352 # gft.trace_configs = [['all', 'all']]
1353 gft.__applyTypicalProcedures()
1354 return gft
1356 def __applyTypicalProcedures(self):
1357 if self.trace_configs is None:
1358 self.createDisplacementTraces()
1359 else:
1360 for src, sen in self.trace_configs:
1361 self.createDisplacementTraces(src, sen)
1362 if self.plot_velocity:
1363 if self.trace_configs is None:
1364 self.createVelocityTraces()
1365 else:
1366 for src, sen in self.trace_configs:
1367 self.createVelocityTraces('{0}|{1}'.format(src, sen))
1368 self.applyFrequencyFilters()
1369 self.getPhaseArrivals()
1371 def __update(self, **kwargs):
1372 for k in kwargs:
1373 temp = kwargs[k]
1374 if temp is not None:
1375 setattr(self, k, temp)
1377 @classmethod
1378 def runStandardCheck(
1379 cls, store_dir, source_depth=None,
1380 lowpass_frequency=None, highpass_frequency=None,
1381 rel_lowpass_frequency=(1. / 110), rel_highpass_frequency=(1. / 16),
1382 distance_min=None, distance_max=None, sensor_count=50,
1383 filter_order=4, pdf_dir=None, plot_velocity=False,
1384 plot_everything=True, output_format='pdf'):
1386 args = locals()
1387 del args['cls']
1388 gft = cls.__createStandardSetups(**args)
1389 gft.createOutputDoc()
1390 gft.__printMessage()
1391 return gft
1393 @classmethod
1394 def runComparissonStandardCheck(
1395 cls, store_dir1, store_dir2, distance_min, distance_max,
1396 together=True, source_depth=None, output_format='pdf',
1397 lowpass_frequency=None, highpass_frequency=None,
1398 rel_lowpass_frequency=(1. / 110), rel_highpass_frequency=(1. / 16),
1399 filter_order=4, pdf_dir=None, plot_velocity=False,
1400 plot_everything=True, sensor_count=50):
1402 args = locals()
1403 del args['cls']
1404 args['change_dists'] = False
1405 gft1 = cls.__createStandardSetups(store_dir1, **args)
1406 if 'source_depth' not in args or args['source_depth'] is None:
1407 args['source_depth'] = gft1.sources[list(gft1.sources.keys())[0]]\
1408 .depth
1410 gft2 = cls.__createStandardSetups(store_dir2, **args)
1412 cls.createComparisonOutputDoc(gft1, gft2, together=together)
1413 gft1.__printMessage()
1414 gft2.__printMessage()
1415 return gft1, gft2
1417 @staticmethod
1418 def createDocumentFromFile(filename, allowed=1, **kwargs):
1419 with open(filename, 'r') as f:
1420 fstr = f.read()
1421 gfts = []
1422 st = None
1423 end = -1
1424 tstr = '--- !{0}.GreensFunctionTest'.format(guts_prefix)
1425 try:
1426 while True:
1427 st = fstr.index(tstr, end + 1)
1428 if st is None:
1429 break
1430 if st > 0 and fstr[st - 1:st] != '\n':
1431 end = st
1432 st = None
1433 continue
1434 end = fstr.index('\n{0}'.format(tstr), st)
1435 gfts.append(fstr[st:end])
1436 except ValueError:
1437 if st is not None:
1438 gfts.append(fstr[st:])
1440 cnt = len(gfts)
1441 if cnt == 0:
1442 raise TypeError('error: GreensFunctionTest: non-valid config'
1443 ' file passed: {0}'.format(filename))
1444 elif cnt > allowed:
1445 raise TypeError('error: GreensFunctionTest: to many configs in'
1446 ' file passed: {0}'.format(filename))
1447 elif cnt < allowed:
1448 raise TypeError('error: GreensFunctionTest: not enough configs in'
1449 ' file passed: {0}'.format(filename))
1451 if cnt == allowed == 1:
1452 gft = load(stream=gfts[0])
1453 gft.__update(**kwargs)
1454 gft.__applyTypicalProcedures()
1455 gft.createOutputDoc()
1456 gft.__printMessage()
1457 return gft
1459 if cnt == allowed == 2:
1460 gft = load(stream=gfts[0])
1461 gft.__update(**kwargs)
1462 gft.__applyTypicalProcedures()
1464 gft2 = load(stream=gfts[1])
1465 gft2.__update(**kwargs)
1466 gft2.__applyTypicalProcedures()
1468 gft.createComparisonOutputDoc(gft, gft2, **kwargs)
1469 gft.__printMessage()
1470 gft2.__printMessage()
1471 return gft, gft2