beach_gmt.py

import gmtpy
import subprocess
import os

wh = 2.0*gmtpy.cm
mw = 0.1*gmtpy.cm

inch = 2.54*gmtpy.cm

beachball_config = dict(
    width=wh,
    height=wh,
    margins=(mw, mw, mw, mw),
    fillcolor=gmtpy.color((255, 0, 0)),
    resolution_dpi=150.,
    pen='1p,black',
    plot_method='sdr',
)


def beachball(mt, filename, **conf_overrides):
    conf = dict(**beachball_config)
    if conf_overrides is not None:
        conf.update(conf_overrides)

    w = conf.pop('width')
    h = conf.pop('height')
    margins = conf.pop('margins')

    indicate_plane = 0
    if 'indicate_plane' in conf:
        indicate_plane = conf.pop('indicate_plane')

    gmtconfig = {}
    for k in conf.keys():
        if k.upper() == k:
            gmtconfig[k] = conf[k]

    gmtconfig['PAPER_MEDIA'] = 'Custom_%ix%i' % (w, h)
    gmtconfig['PS_MITER_LIMIT'] = '180'

    gmt = gmtpy.GMT(config=gmtconfig)

    layout = gmt.default_layout()
    layout.set_fixed_margins(*margins)

    widget = layout.get_widget()

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

    m = mt.m_up_south_east() / mt.scalar_moment() * 4.
    mrr = m[0, 0]
    mtt = m[1, 1]
    mff = m[2, 2]
    mrt = m[0, 1]
    mrf = m[0, 2]
    mtf = m[1, 2]

    if conf['plot_method'] == 'sdr':
        kwargs = dict(R=(-1., 1., -1., 1.),
                      S='a%gp' % ((w-margins[0]-margins[1])),
                      in_rows=[[0., 0., 1.,
                                strike, dip, rake,
                                5.0, 0., 0., '']])

    elif conf['plot_method'] == 'zerotrace':
        kwargs = dict(R=(-1., 1., -1., 1.),
                      S='m%gp' % ((w-margins[0]-margins[1])),
                      in_rows=[[0., 0., 1.,
                                mrr, mtt, mff, mrt, mrf, mtf,
                                23., 0., 0., '']])
    else:
        assert False

    if indicate_plane == 0:
        gmt.psmeca(N=True, L=conf['pen'], G=conf['fillcolor'], *widget.JXY(),
                   **kwargs)
    else:
        gmt.psmeca(N=True, L='1p,black,.', W=conf['pen'], G=conf['fillcolor'],
                   *widget.JXY(), **kwargs)
        gmt.psmeca(N=True, T='%i/2p,%s' % (indicate_plane,
                                           gmtpy.color('black')),
                   *widget.JXY(), **kwargs)

    if filename.endswith('.png'):
        gmt.save(filename+'.pdf')
        subprocess.check_call(
            ['pdftocairo', '-singlefile', '-r', str(conf['resolution_dpi']),
             '-png', filename+'.pdf', filename[:-4]])
        os.remove(filename+'.pdf')

    return filename