Plotting functions

Generating topographic maps with automap

The pyrocko.plot.automap module provides a painless and clean interface for the Generic Mapping Tool (GMT) [1].

Classes covered in these examples:

• pyrocko.plot.automap.Map

For details on our approach in calling GMT from Python:

• GmtPy - A Python interface to GMT 4/5

Note

To retain PDF transparency in gmtpy use save(psconvert=True).

Topographic map of Dead Sea basin

This example demonstrates how to create a map of the Dead Sea area with largest cities, topography and gives a hint on how to access genuine GMT methods.

Download automap_example.py

Station file used in the example: stations_deadsea.pf

from pyrocko.plot.automap import Map
from pyrocko.example import get_example_data
from pyrocko import model
from pyrocko import moment_tensor as pmt
from pyrocko.plot import gmtpy
gmtpy.check_have_gmt()

# Download example data
get_example_data('stations_deadsea.pf')
get_example_data('deadsea_events_1976-2017.txt')

# Generate the basic map
m = Map(
    lat=31.5,
    lon=35.5,
    radius=150000.,
    width=30., height=30.,
    show_grid=False,
    show_topo=True,
    color_dry=(238, 236, 230),
    topo_cpt_wet='light_sea_uniform',
    topo_cpt_dry='light_land_uniform',
    illuminate=True,
    illuminate_factor_ocean=0.15,
    show_rivers=False,
    show_plates=True)

# Draw some larger cities covered by the map area
m.draw_cities()

# Generate with latitute, longitude and labels of the stations
stations = model.load_stations('stations_deadsea.pf')
lats = [s.lat for s in stations]
lons = [s.lon for s in stations]
labels = ['.'.join(s.nsl()) for s in stations]

# Stations as black triangles. Genuine GMT commands can be parsed by the maps'
# gmt attribute. Last argument of the psxy function call pipes the maps'
# pojection system.
m.gmt.psxy(in_columns=(lons, lats), S='t20p', G='black', *m.jxyr)

# Station labels
for i in range(len(stations)):
    m.add_label(lats[i], lons[i], labels[i])


# Load events from catalog file (generated using catalog.GlobalCMT()
# download from www.globalcmt.org)
# If no moment tensor is provided in the catalogue, the event is plotted
# as a red circle. Symbol size relative to magnitude.

events = model.load_events('deadsea_events_1976-2017.txt')
beachball_symbol = 'd'
factor_symbl_size = 5.0
for ev in events:
    mag = ev.magnitude
    if ev.moment_tensor is None:
        ev_symb = 'c'+str(mag*factor_symbl_size)+'p'
        m.gmt.psxy(
            in_rows=[[ev.lon, ev.lat]],
            S=ev_symb,
            G=gmtpy.color('scarletred2'),
            W='1p,black',
            *m.jxyr)
    else:
        devi = ev.moment_tensor.deviatoric()
        beachball_size = mag*factor_symbl_size
        mt = devi.m_up_south_east()
        mt = mt / ev.moment_tensor.scalar_moment() \
            * pmt.magnitude_to_moment(5.0)
        m6 = pmt.to6(mt)
        data = (ev.lon, ev.lat, 10) + tuple(m6) + (1, 0, 0)

        if m.gmt.is_gmt5():
            kwargs = dict(
                M=True,
                S='%s%g' % (beachball_symbol[0], (beachball_size) / gmtpy.cm))
        else:
            kwargs = dict(
                S='%s%g' % (beachball_symbol[0],
                            (beachball_size)*2 / gmtpy.cm))

        m.gmt.psmeca(
            in_rows=[data],
            G=gmtpy.color('chocolate1'),
            E='white',
            W='1p,%s' % gmtpy.color('chocolate3'),
            *m.jxyr,
            **kwargs)

m.save('automap_deadsea.png')

Footnotes

[1]	Wessel, P., W. H. F. Smith, R. Scharroo, J. F. Luis, and F. Wobbe, Generic Mapping Tools: Improved version released, EOS Trans. AGU, 94, 409-410, 2013.

Plotting beachballs (focal mechanisms)

Classes covered in these examples:

• pyrocko.plot.beachball (visual representation of a focal mechanism)
• pyrocko.moment_tensor (a 3x3 matrix representation of an earthquake source)
• pyrocko.gf.seismosizer.DCSource (a representation of a double couple source object),
• pyrocko.gf.seismosizer.RectangularExplosionSource (a representation of a rectangular explosion source),
• pyrocko.gf.seismosizer.CLVDSource (a representation of a compensated linear vector dipole source object)
• pyrocko.gf.seismosizer.DoubleDCSource (a representation of a double double-couple source object).

Beachballs from moment tensors

Here we create random moment tensors and plot their beachballs.

Download beachball_example01.py

import random
import logging
from matplotlib import pyplot as plt
from pyrocko import moment_tensor as pmt
from pyrocko import util
from pyrocko.plot import beachball

''' Beachball Copacabana '''

logger = logging.getLogger('pyrocko.examples.beachball_example01')

util.setup_logging()

fig = plt.figure(figsize=(10., 4.))
fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
axes = fig.add_subplot(1, 1, 1)

for i in range(200):

    # create random moment tensor
    mt = pmt.MomentTensor.random_mt()

    try:
        # create beachball from moment tensor
        beachball.plot_beachball_mpl(
            mt, axes,
            # type of beachball: deviatoric, full or double couple (dc)
            beachball_type='full',
            size=random.random()*120.,
            position=(random.random()*10., random.random()*10.),
            alpha=random.random(),
            linewidth=1.0)

    except beachball.BeachballError as e:
        logger.error('%s for MT:\n%s' % (e, mt))

axes.set_xlim(0., 10.)
axes.set_ylim(0., 10.)
axes.set_axis_off()
fig.savefig('beachball-example01.pdf')

plt.show()

An artistic display of focal mechanisms drawn by classes pyrocko.plot.beachball and pyrocko.moment_tensor.

This example shows how to plot a full, a deviatoric and a double-couple beachball for a moment tensor.

Download beachball_example03.py

from matplotlib import pyplot as plt
from pyrocko import moment_tensor as pmt
from pyrocko import plot
from pyrocko.plot import beachball


fig = plt.figure(figsize=(4., 2.))
fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
axes = fig.add_subplot(1, 1, 1)
axes.set_xlim(0., 4.)
axes.set_ylim(0., 2.)
axes.set_axis_off()

for i, beachball_type in enumerate(['full', 'deviatoric', 'dc']):
    beachball.plot_beachball_mpl(
            pmt.as_mt((124654616., 370943136., -6965434.0,
                       553316224., -307467264., 84703760.0)),
            axes,
            beachball_type=beachball_type,
            size=60.,
            position=(i+1, 1),
            color_t=plot.mpl_color('scarletred2'),
            linewidth=1.0)

fig.savefig('beachball-example03.pdf')
plt.show()

The three types of beachballs that can be plotted through pyrocko.

Beachballs from source objects

This example shows how to add beachballs of various sizes to the corners of a plot by obtaining the moment tensor from four different source object types: pyrocko.gf.seismosizer.DCSource (upper left), pyrocko.gf.seismosizer.RectangularExplosionSource (upper right), pyrocko.gf.seismosizer.CLVDSource (lower left) and pyrocko.gf.seismosizer.DoubleDCSource (lower right).

Creating the beachball this ways allows for finer control over their location based on their size (in display units) which allows for a round beachball even if the axis are not 1:1.

Download beachball_example02.py

from matplotlib import transforms, pyplot as plt
from pyrocko import gf
from pyrocko.plot import beachball


# create source object
source1 = gf.DCSource(depth=35e3, strike=0., dip=90., rake=0.)

# set size of beachball
sz = 20.
# set beachball offset in points (one point from each axis)
szpt = (sz / 2.) / 72. + 1. / 72.

fig = plt.figure(figsize=(10., 4.))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlim(0, 10)
ax.set_ylim(0, 10)

# get the bounding point (left-top)
x0 = ax.get_xlim()[0]
y1 = ax.get_ylim()[1]

# create a translation matrix, based on the final figure size and
# beachball location
move_trans = transforms.ScaledTranslation(szpt, -szpt, fig.dpi_scale_trans)

# get the inverse matrix for the axis where the beachball will be plotted
inv_trans = ax.transData.inverted()

# set the bouding point relative to the plotted axis of the beachball
x0, y1 = inv_trans.transform(move_trans.transform(
    ax.transData.transform((x0, y1))))

# plot beachball
beachball.plot_beachball_mpl(source1.pyrocko_moment_tensor(), ax,
                             beachball_type='full', size=sz,
                             position=(x0, y1), linewidth=1.)


# create source object
source2 = gf.RectangularExplosionSource(depth=35e3, strike=0., dip=90.)

# set size of beachball
sz = 30.
# set beachball offset in points (one point from each axis)
szpt = (sz / 2.) / 72. + 1. / 72.

# get the bounding point (right-upper)
x1 = ax.get_xlim()[1]
y1 = ax.get_ylim()[1]

# create a translation matrix, based on the final figure size and
# beachball location
move_trans = transforms.ScaledTranslation(-szpt, -szpt, fig.dpi_scale_trans)

# get the inverse matrix for the axis where the beachball will be plotted
inv_trans = ax.transData.inverted()

# set the bouding point relative to the plotted axis of the beachball
x1, y1 = inv_trans.transform(move_trans.transform(
    ax.transData.transform((x1, y1))))

# plot beachball
beachball.plot_beachball_mpl(source2.pyrocko_moment_tensor(), ax,
                             beachball_type='full', size=sz,
                             position=(x1, y1), linewidth=1.)


# create source object
source3 = gf.CLVDSource(moment=1.0, azimuth=30., dip=30.)

# set size of beachball
sz = 40.
# set beachball offset in points (one point from each axis)
szpt = (sz / 2.) / 72. + 1. / 72.

# get the bounding point (left-bottom)
x0 = ax.get_xlim()[0]
y0 = ax.get_ylim()[0]

# create a translation matrix, based on the final figure size and
# beachball location
move_trans = transforms.ScaledTranslation(szpt, szpt, fig.dpi_scale_trans)

# get the inverse matrix for the axis where the beachball will be plotted
inv_trans = ax.transData.inverted()

# set the bouding point relative to the plotted axis of the beachball
x0, y0 = inv_trans.transform(move_trans.transform(
    ax.transData.transform((x0, y0))))

# plot beachball
beachball.plot_beachball_mpl(source3.pyrocko_moment_tensor(), ax,
                             beachball_type='full', size=sz,
                             position=(x0, y0), linewidth=1.)

# create source object
source4 = gf.DoubleDCSource(depth=35e3, strike1=0., dip1=90., rake1=0.,
                            strike2=45., dip2=74., rake2=0.)

# set size of beachball
sz = 50.
# set beachball offset in points (one point from each axis)
szpt = (sz / 2.) / 72. + 1. / 72.

# get the bounding point (right-bottom)
x1 = ax.get_xlim()[1]
y0 = ax.get_ylim()[0]

# create a translation matrix, based on the final figure size and
# beachball location
move_trans = transforms.ScaledTranslation(-szpt, szpt, fig.dpi_scale_trans)

# get the inverse matrix for the axis where the beachball will be plotted
inv_trans = ax.transData.inverted()

# set the bouding point relative to the plotted axis of the beachball
x1, y0 = inv_trans.transform(move_trans.transform(
    ax.transData.transform((x1, y0))))

# plot beachball
beachball.plot_beachball_mpl(source4.pyrocko_moment_tensor(), ax,
                             beachball_type='full', size=sz,
                             position=(x1, y0), linewidth=1.)

fig.savefig('beachball-example02.pdf')
plt.show()

Four different source object types plotted with different beachball sizes.

Fuzzy beachballs with uncertainty

If we want to express moment tensor uncertainties we can plot fuzzy beachballs from an ensemble of many solutions.

This example will generate random solution around a best moment tensor (red lines). The perturbed solutions are the uncertainty which can be illustrated in a fuzzy beachball.

Download beachball_example05.py

from pyrocko.plot import beachball
import pyrocko.moment_tensor as mtm

import numpy as num
from matplotlib import pyplot as plt


fig = plt.figure(figsize=(4., 4.))
fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
axes = fig.add_subplot(1, 1, 1)

# Number of available solutions
n_balls = 1000

# Best solution
strike = 135.
dip = 65.
rake = 15.

best_mt = mtm.MomentTensor.from_values((strike, dip, rake))


def get_random_uniform(lower, upper, dimension=1):
    ''' Help function to pertubate the best solution '''
    return (upper - lower) * num.random.rand(dimension) + lower


mts = []
for i in range(n_balls):
    strike_dev = get_random_uniform(-15., 15.)
    mts.append(mtm.MomentTensor.from_values(
        (strike + strike_dev, dip, rake)))

plot_kwargs = {
    'beachball_type': 'full',
    'size': 8,
    'position': (5, 5),
    'color_t': 'black',
    'edgecolor': 'black'
    }

beachball.plot_fuzzy_beachball_mpl_pixmap(mts, axes, best_mt, **plot_kwargs)

# Decorate the axes
axes.set_xlim(0., 10.)
axes.set_ylim(0., 10.)
axes.set_axis_off()

plt.show()

Fuzzy beachball illustrating the solutions uncertainty.

Beachballs views for cross-sections:

It is useful to show beachballs from other view angles, as in cross-sections. For that, we can define a view for all beachball plotting functions as shown here:

Download beachball_example06.py

#!/usr/bin/env python3

from matplotlib import pyplot as plt
from pyrocko import moment_tensor as pmt
from pyrocko.plot import beachball

mt = pmt.as_mt([0.424, -0.47, 0.33, 0.711, -0.09, 0.16])
axes = plt.gca()

beachball.plot_beachball_mpl(
    mt, axes,
    size=50.,
    position=(0., 0.),
    view='top')

beachball.plot_beachball_mpl(
    mt, axes,
    size=50.,
    position=(0, -1.),
    view='south')

beachball.plot_beachball_mpl(
    mt, axes,
    size=50.,
    position=(-1, 0.),
    view='east')

beachball.plot_beachball_mpl(
    mt, axes,
    size=50.,
    position=(0, 1.),
    view='north')

beachball.plot_beachball_mpl(
    mt, axes,
    size=50.,
    position=(1, 0.),
    view='west')

axes.set_xlim(-2., 2.)
axes.set_ylim(-2., 2.)
plt.show()

Add station symbols to focal sphere diagram

This example shows how to add station symbols at the positions where P wave rays pierce the focal sphere.

The function to plot focal spheres (pyrocko.plot.beachball.plot_beachball_mpl()) uses the function pyrocko.plot.beachball.project() in the final projection from 3D to 2D coordinates. Here we use this function to place additional symbols on the plot. The take-off angles needed can be computed with some help of the pyrocko.cake module. Azimuth and distance computations are done with functions from pyrocko.orthodrome.

Download beachball_example04.py

import numpy as num
from matplotlib import pyplot as plt
from pyrocko import moment_tensor as pmt, cake, orthodrome
from pyrocko.plot import beachball

km = 1000.

# source position and mechanism
slat, slon, sdepth = 0., 0., 10.*km
mt = pmt.MomentTensor.random_dc()

# receiver positions
rdepth = 0.0
rlatlons = [(50., 10.), (60., -50.), (-30., 60.)]

# earth model and phase for takeoff angle computations
mod = cake.load_model('ak135-f-continental.m')
phases = cake.PhaseDef.classic('P')

# setup figure with aspect=1.0/1.0, ranges=[-1.1, 1.1]
fig = plt.figure(figsize=(2., 2.))  # size in inch
fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
axes = fig.add_subplot(1, 1, 1, aspect=1.0)
axes.set_axis_off()
axes.set_xlim(-1.1, 1.1)
axes.set_ylim(-1.1, 1.1)

projection = 'lambert'

beachball.plot_beachball_mpl(
    mt, axes,
    position=(0., 0.),
    size=2.0,
    color_t=(0.7, 0.4, 0.4),
    projection=projection,
    size_units='data')

for rlat, rlon in rlatlons:
    distance = orthodrome.distance_accurate50m(slat, slon, rlat, rlon)
    rays = mod.arrivals(
        phases=cake.PhaseDef('P'),
        zstart=sdepth, zstop=rdepth, distances=[distance*cake.m2d])

    if not rays:
        continue

    takeoff = rays[0].takeoff_angle()
    azi = orthodrome.azimuth(slat, slon, rlat, rlon)

    # to spherical coordinates, r, theta, phi in radians
    rtp = num.array([[1., num.deg2rad(takeoff), num.deg2rad(90.-azi)]])

    # to 3D coordinates (x, y, z)
    points = beachball.numpy_rtp2xyz(rtp)

    # project to 2D with same projection as used in beachball
    x, y = beachball.project(points, projection=projection).T

    axes.plot(x, y, '+', ms=10., mew=2.0, mec='black', mfc='none')

fig.savefig('beachball-example04.png')

Focal sphere diagram with markers at positions of P wave ray piercing points.

Hudson’s source type plot

Hudson’s source type plot [Hudson, 1989] is a way to visually represent the widely used “standard” decomposition of a moment tensor into its isotropic, its compensated linear vector dipole (CLVD), and its double-couple (DC) components.

The function pyrocko.plot.hudson.project() may be used to get the (u,v) coordinates for a given (full) moment tensor used for positioning the symbol in the plot. The function pyrocko.plot.hudson.draw_axes() can be used to conveniently draw the axes and annotations. Note, that we follow the original convention introduced by Hudson, to place the negative CLVD on the right hand side.

Download hudson_diagram.py

from __future__ import print_function
import sys
from matplotlib import pyplot as plt
from pyrocko.plot import hudson, beachball, mpl_init, mpl_color
from pyrocko import moment_tensor as pmt

# a bunch of random MTs
moment_tensors = [pmt.random_mt() for _ in range(200)]

# setup plot layout
fontsize = 10.
markersize = fontsize
mpl_init(fontsize=fontsize)
width = 7.
figsize = (width, width / (4. / 3.))
fig = plt.figure(figsize=figsize)
axes = fig.add_subplot(1, 1, 1)
fig.subplots_adjust(left=0.03, right=0.97, bottom=0.03, top=0.97)

# draw focal sphere diagrams for the random MTs
for mt in moment_tensors:
    u, v = hudson.project(mt)
    try:
        beachball.plot_beachball_mpl(
            mt, axes,
            beachball_type='full',
            position=(u, v),
            size=markersize,
            color_t=mpl_color('skyblue3'),
            color_p=mpl_color('skyblue1'),
            alpha=1.0,  # < 1 for transparency
            zorder=1,
            linewidth=0.25)

    except beachball.BeachballError as e:
        print(str(e), file=sys.stderr)

# draw the axes and annotations of the hudson plot
hudson.draw_axes(axes)

fig.savefig('hudson_diagram.png', dpi=150)
# plt.show()

Hudson’s source type plot for 200 random moment tensors.

Source radiation plot

The directivity and radiation characteristics of any point or finite Source model can be illustrated with plot_directivity().

Radiation pattern effects

The following educational example illustrates radiation pattern effects from a point source in a homogeneous full space. Analytical Green’s functions for a homogeneous full space are computed within the example code by use of the ahfullgreen backend of Fomosto.

Download plot_directivity.py

import os
import shutil
from pyrocko.plot.directivity import plot_directivity
from pyrocko.gf import LocalEngine, DCSource, Store
from pyrocko.fomosto import ahfullgreen

km = 1e3


def make_homogeneous_gf_store(
        path, store_id, source_depth, receiver_depth, distance):

    if os.path.exists(path):
        shutil.rmtree(path)

    ahfullgreen.init(path, None, config_params=dict(
        id=store_id,
        sample_rate=20.,
        receiver_depth=receiver_depth,
        source_depth_min=source_depth,
        source_depth_max=source_depth,
        distance_min=distance,
        distance_max=distance))

    store = Store(path)
    store.make_ttt()
    ahfullgreen.build(path)


store_id = 'gf_homogeneous_radpat'
store_path = os.path.join('.', store_id)
distance = 10*km
receiver_depth = 0.0

source = DCSource(
    depth=0.,
    strike=0.,
    dip=90.,
    rake=0.)

make_homogeneous_gf_store(
    store_path, store_id, source.depth, receiver_depth, distance)

engine = LocalEngine(store_dirs=[store_path])

# import matplotlib.pyplot as plt
# fig = plt.figure(figsize=(7,5))
# axes = fig.add_subplot(111, polar=True)

resp = plot_directivity(
    engine, source, store_id,
    # axes=axes,
    distance=distance,
    dazi=5.,
    component='R',
    target_depth=receiver_depth,
    plot_mt='full',
    show_phases=True,
    fmin=None,
    fmax=1.0,
    phases={
        'P': '{stored:anyP}-50%',
        'S': '{stored:anyS}+50%'
    },
    interpolation='nearest_neighbor',
    quantity='velocity',
    envelope=False,
    hillshade=False)

# fig.savefig('radiation_pattern.png')

Radial component radiation pattern for a double-couple point source in a homogeneous full space observed in the plane of the source at a distance of 10 km. Note that the S waves seen in this example are pure near-field effects. They get less pronounced when going to higher frequencies.

Directivity effects

Synthetic seismic traces (R, T or Z) are forward-modelled at a defined radius, covering the full or partial azimuthal range and projected on a polar plot. Difference in the amplitude are enhanced by hillshading the data.

Download plot_directivity.py

import os
from pyrocko.plot.directivity import plot_directivity
from pyrocko.gf import LocalEngine, RectangularSource, ws

km = 1e3
# The store we are going extract data from:
store_id = 'iceland_reg_v2'

# First, download a Greens Functions store. If you already have one that you
# would like to use, you can skip this step and point the *store_superdirs* in
# the next step to that directory.

if not os.path.exists(store_id):
    ws.download_gf_store(site='kinherd', store_id=store_id)

# We need a pyrocko.gf.Engine object which provides us with the traces
# extracted from the store.
engine = LocalEngine(store_superdirs=['.'])

# Create a RectangularSource with uniform fit.
rect_source = RectangularSource(
    depth=1.6*km,
    strike=240.,
    dip=76.6,
    rake=-.4,
    anchor='top',
    nucleation_x=-.57,
    nucleation_y=-.59,
    velocity=2070.,
    length=27*km,
    width=9.4*km,
    slip=1.4)


# import matplotlib.pyplot as plt
# fig = plt.figure(figsize=(7,5))
# axes = fig.add_subplot(111, polar=True)

resp = plot_directivity(
    engine, rect_source, store_id,
    # axes=axes,
    distance=300*km,
    dazi=5.,
    component='R',
    plot_mt='full',
    show_phases=True,
    phases={
        'First': 'first{stored:begin}-10%',
        'Last': 'last{stored:end}+20'
    },
    quantity='displacement',
    envelope=False)

# fig.savefig('directivity_rectangular.png')

Source radiation pattern at 300 km distance of the Mw 6.8 2020 Elazig-Sevrice earthquake. The dominantly unilateral strike-slip rupture is reconstructed by a finite RectangularSource model.

Here we see the envelope of the synthetic seismic traces, emphasizing the directivity effects of the source (envelope=True). Same source model: Mw 6.8 2020 Elazig-Sevrice earthquake.