bb_beauty_contest.py

from subprocess import check_call, check_output
import os
import random
import shutil
import traceback
import tempfile
from pyrocko import moment_tensor as pmt
from obspy.imaging import beachball as opb
import gmtpy
import mopad
import beach_gmt
import numpy as num

inch2cm = 2.54


def tempfn(*args):
    return tempfile.mkstemp(*args)[1]


def m6(m):
    return m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2]


def plot_mopad(mt, fn_out):

    mt = mopad.MomentTensor(M=mt.m())
    beach = mopad.BeachBall(mt)
    fn = tempfn('.png')
    beach.save_BB(dict(
        plot_outfile=fn,
        plot_outfile_format='png',
        plot_size=2.,
        plot_dpi=100,
        plot_projection='lambert',
        plot_nodalline_width=5.,
        plot_outerline_width=5.))

    check_call(['convert', '-flatten', fn, fn_out])
    os.remove(fn)


def plot_gmt(mt, fn_out, **kwargs):
    size = 2*inch2cm*gmtpy.cm
    beach_gmt.beachball(
        mt, fn_out,
        resolution_dpi=100.,
        width=size, height=size,
        margins=[size/40.]*4,
        pen='5p,black', **kwargs)


def plot_obspy(mt, fn_out):
    m = mt.m_up_south_east()
    err = None
    fn = tempfn('.png')
    try:
        opb.Beachball(m6(m), linewidth=5, width=200.,
                      facecolor='r',
                      outfile=fn)
        check_call(['convert', '-flatten', fn, fn_out])

    except Exception:
        err = traceback.format_exc()
        check_call(['convert', '-size', '200x200', 'xc:white', fn_out])

    finally:
        os.remove(fn)

    return err


def fuzz(mi, ma, critical=[]):

    random_amp = random.choice([0.5, 1e-3, 1e-6, 1e-12])

    r = random.random() * random_amp
    m = (mi+ma)*0.5

    candidates = [mi, ma, mi + r, ma - r, m, m - r, m + r]

    for c in critical:
        for x in (c, c-r, c+r):
            if mi <= x <= ma:
                candidates.append(x)

    return random.choice(candidates)


def random_mt(kind):
    if kind == 'dc':
        mt = pmt.MomentTensor.random_dc()
        strike, dip, rake = mt.both_strike_dip_rake()[0]

    elif kind == 'zerotrace':
        m = pmt.symmat6(*num.random.random(6)*2.-1.)
        t = float(m.trace()/3.)
        m -= num.diag([t, t, t])
        mt = pmt.MomentTensor(m=m)

    elif kind == 'critical-dc':
        strike = fuzz(-180., 180., critical=[-90., 90.])
        dip = fuzz(0., 90.)
        rake = fuzz(-180., 180., critical=[-90., 90.])
        mt = pmt.MomentTensor(strike=strike, dip=dip, rake=rake)

        for sdr in mt.both_strike_dip_rake():
            mt2 = pmt.MomentTensor(strike=sdr[0], dip=sdr[1], rake=sdr[2])
            assert num.max(num.abs(mt.m() - mt2.m())) < 1e-7

    return mt


def valform(*vals):
    fmt = ', '.join(['%.15f'] * len(vals))
    return fmt % tuple(vals)

testkinds = 'critical-dc'
ifail = 0
o = open('comparison.html', 'w')
for i in xrange(50):

    mt = random_mt(testkinds)
    strike, dip, rake = mt.both_strike_dip_rake()[0]

    fn_beach = {}
    fn = fn_beach['mopad'] = tempfn('.png')
    plot_mopad(mt, fn)

    fn = fn_beach['gmt'] = tempfn('.png')
    plot_gmt(mt, fn, plot_method='zerotrace')

    fn = fn_beach['obspy'] = tempfn('.png')
    err = plot_obspy(mt, fn)

    diff = {}
    fns_comp = []
    for a, b in [('mopad', 'gmt'), ('mopad', 'obspy')]:
        #compfn = 'compare-%s-%s.png' % (a, b)
        compfn = tempfn('.png')

        check_call(['convert', fn_beach[a], fn_beach[b], '-compose',
                    'difference', '-composite', '-colorspace', 'gray', compfn])
        diff[a, b] = float(check_output(['convert', compfn, '-format',
                                         '%[fx:mean*100]', 'info:']))
        fns_comp.append(compfn)

    fn_result = tempfn('.png')
    fns_in = [fn_beach[x] for x in ['mopad', 'gmt', 'obspy']]
    check_call(['convert'] + fns_in +
               fns_comp + ['+append', fn_result])

    for fn in fns_comp + fns_in:
        os.remove(fn)

    if any(v > 4.0 for v in diff.values()):
        ifail += 1
        shutil.copy(fn_result, 'fail-%i.png' % ifail)
        for ((a, b), v) in diff.iteritems():
            print >>o, 'difference %s vs %s: %g%%<br />' % (a, b, v)

        print >>o, 'strike, dip, rake: %s<br />' % valform(strike, dip, rake)
        print >>o, 'mnn, mee, mdd, mne, mnd, med: %s<br />' % \
            valform(*m6(mt.m()))

        if err is not None:
            print >>o, 'obspy failure: <pre>%s</pre><br />' % err

        print >>o, '<img src="fail-%i.png"> <br />' % ifail
        print >>o, '(mopad, gmt, obspy, |gmt - mopad|, |obpsy - mopad|) <br />'
        print >>o, '<hr />'

    os.remove(fn_result)