WIP: Add tutorial showing how to plot focal mechanisms (beachballs)

examples/tutorials/advanced/focal_mechanisms.py

@@ -0,0 +1,370 @@
+"""
+Plotting focal mechanisms
+=========================
+
+Focal mechanisms can be plotted as beachballs with the :meth:`pygmt.Figure.meca`
+method.
+
+The focal mechanism data or parameters can be provided as different input
+types: external file, 1-D or 2-D ``numpy.array``, dictionary, ``pandas.Dataframe``.
+Different conventions to define the focal mechanism are supported: Aki and
+Richards (``"aki"``), global CMT (``"gcmt"``), moment tensor (``"mt"``), partial
+focal mechanism (``"partial"``), principal axis (``"principal_axis"``). Please
+refer to the documentation of :meth:`pygmt.Figure.meca` regarding how to set up
+the input data in respect to the chosen input type and convention (i.e the
+expected column order, keys, or column names).
+
+This tutorial focus on how to adjust the display of the beachballs:
+
+- Filling the quadrants
+- Plotting the P and T axes
+- Plotting the components of a seismic moment tensor
+- Adjusting the outlines
+- Highlighting the nodal planes
+- Adding offset from the event location
+- Adding a label
+- Using size-coding and color-coding
+"""
+
+# %%
+import pandas as pd
+import pygmt
+
+# Set up arguments for basemap
+size = 5
+projection = "X10c/4c"
+frame = ["af", "+ggray90"]
+
+
+# %%
+# Set up the focal mechanism data
+# -------------------------------
+#
+# Store focal mechanism parameters for one event in a dictionary following the
+
+# moment tensor convention
+mt_dict_single = mt_virginia = {
+    "mrr": 4.71,
+    "mtt": 0.0381,
+    "mff": -4.74,
+    "mrt": 0.399,
+    "mrf": -0.805,
+    "mtf": -1.23,
+    "exponent": 24,
+}
+# Aki and Richards convention
+aki_dict_single = {"strike": 318, "dip": 89, "rake": -179, "magnitude": 7.75}
+
+
+# %%
+# Plotting a single beachball
+# ---------------------------
+#
+# Required parameters are ``spec``, ``scale``, ``longitude`` and ``latitude``
+# (event location) as well as ``convention``. For the input types dictionary and
+# ``pandas.Dataframe``, ``convention`` is not required.
+
+fig = pygmt.Figure()
+fig.basemap(region=[-size, size] * 2, projection=projection, frame=frame)
+
+fig.meca(spec=mt_dict_single, scale="1c", longitude=0, latitude=0)
+
+fig.show()
+
+
+# %%
+# Filling the quadrants
+# ---------------------
+#
+# Use the parameters ``compressionfill`` and ``extensionfill`` to fill the
+# quadrants with different colors or patterns. Regarding pattern see the
+# gallery example :doc:`Bit and hachure patterns </gallery/symbols/patterns>`
+# and the Technical Reference :doc:`Bit and hachure patterns </techref/patterns>`.
+
+fig = pygmt.Figure()
+fig.basemap(region=[-size, size] * 2, projection=projection, frame=frame)
+
+fig.meca(
+    spec=mt_dict_single,
+    scale="1c",
+    longitude=-2,
+    latitude=0,
+    compressionfill="darkorange",
+    extensionfill="cornsilk",
+)
+
+fig.meca(
+    spec=mt_dict_single,
+    scale="1c",
+    longitude=2,
+    latitude=0,
+    compressionfill="p8",
+    extensionfill="p31",
+    outline=True,
+)
+
+fig.show()
+
+
+# %%
+# Plotting the P and T axes
+# -------------------------
+#
+# Wait for PR #3526
+
+fig = pygmt.Figure()
+fig.basemap(region=[-size, size] * 2, projection=projection, frame=frame)
+
+fig.meca(
+    spec=mt_dict_single,
+    scale="1c",
+    longitude=-2,
+    latitude=0,
+    Fa=True,
+)
+
+fig.meca(
+    spec=mt_dict_single,
+    scale="1c",
+    longitude=2,
+    latitude=0,
+    Fa="0.2c/cd",  # Compute and plot P and T axes with symbols, Adjust size and symbols
+    Fe="darkorange",  # Adjust fill of T axis symbol
+    Fg="gray30",  # Adjust fill of P axis symbol
+    Ft="0.8p,cornsilk",  # Adjust outline of T axis symbol
+    Fp="0.8p,gray60",  # Adjust outline of P axis symbol
+)
+
+fig.show()
+
+
+# %%
+# Plotting the components of a seismic moment tensor
+# --------------------------------------------------
+#
+# Use the ``component`` parameter to plot the components of a seismic moment tensor.
+
+fig = pygmt.Figure()
+fig.basemap(region=[-size, size] * 2, projection=projection, frame=frame)
+
+fig.meca(
+    spec=mt_dict_single,
+    scale="1c",
+    longitude=-3,
+    latitude=0,
+    component="full",  # full seismic moment tensor
+)
+fig.meca(
+    spec=mt_dict_single,
+    scale="1c",
+    longitude=0,
+    latitude=0,
+    component="dc",  # closest double couple
+)
+fig.meca(
+    spec=mt_dict_single,
+    scale="1c",
+    longitude=2,
+    latitude=0,
+    component="deviatoric",  # deviatoric part
+)
+
+fig.show()
+
+
+# %%
+# Adjusting the outlines
+# ----------------------
+#
+# Use the parameters ``pen`` and ``outline`` for adjusting the circumference of
+# the beachball or all lines (circumference of the beachball and both nodal planes).
+
+fig = pygmt.Figure()
+fig.basemap(region=[-size, size] * 2, projection=projection, frame=frame)
+
+fig.meca(
+    spec=aki_dict_single,
+    scale="1c",
+    longitude=-2,
+    latitude=0,
+    # Use a 1-point thick, darkorange and solid line
+    pen="1p,darkorange,solid",
+)
+
+fig.meca(
+    spec=aki_dict_single,
+    scale="1c",
+    longitude=2,
+    latitude=0,
+    outline="1p,darkorange,solid",
+)
+
+fig.show()
+
+
+# %%
+# Highlighting the nodal planes
+# -----------------------------
+#
+# Use the parameter ``nodal``, whereby ``"0"`` refers to both, ``"1"`` to the
+# first, and ``"2"`` to the second nodal plane(s). Only the circumference and the
+# specified nodal plane(s) are plotted, i.e. the quadrants remain unfilled
+# (transparent). If needed, make usage of the stacking concept of (Py)GMT and use
+# ``nodal`` with the ``outline`` or / and ``pen``  parameters in combination.
+
+fig = pygmt.Figure()
+fig.basemap(region=[-size, size] * 2, projection=projection, frame=frame)
+
+fig.meca(
+    spec=aki_dict_single,
+    scale="1c",
+    longitude=-2,
+    latitude=0,
+    nodal="0/1p,black,solid",
+)
+
+fig.meca(
+    spec=aki_dict_single,
+    scale="1c",
+    longitude=2,
+    latitude=0,
+    compressionfill="lightorange",
+    outline="0.5p,black,solid",
+)
+fig.meca(
+    spec=aki_dict_single,
+    scale="1c",
+    longitude=2,
+    latitude=0,
+    nodal="1/1p,darkorange,solid",
+)
+fig.meca(
+    spec=aki_dict_single,
+    scale="1c",
+    longitude=2,
+    latitude=0,
+    compressionfill="white@100",
+    extensionfill="white@100",
+    pen="1p,gray30,solid",
+)
+fig.show()
+
+
+# %%
+# Adding offset from event location
+# ---------------------------------
+#
+# Specify the optional parameters ``plot_longitude`` and ``plot_latitude``.
+# Additionally the parameter ``offset`` has to be set. Besides just drawing a
+# line between the beachball and the event location, a small circle can be
+# plotted at the event location by appending **+s** and the descired circle
+# diameter. The connecting line as well as the outline of the circle are
+# plotted with the setting of pen, or can be adjusted separately. The fill of
+# the small circle corresponds to the fill of the compressive quadrantes.
+
+fig = pygmt.Figure()
+fig.basemap(region=[-size, size] * 2, projection=projection, frame=frame)
+
+fig.meca(
+    spec=aki_dict_single,
+    scale="1c",
+    longitude=-1,
+    latitude=0,
+    plot_longitude=-3,
+    plot_latitude=2,
+    offset=True,
+)
+
+fig.meca(
+    spec=aki_dict_single,
+    scale="1c",
+    longitude=3,
+    latitude=0,
+    plot_longitude=1,
+    plot_latitude=2,
+    offset="+p1p,darkorange+s0.25c",
+    compressionfill="lightorange",
+)
+
+fig.show()
+
+
+# %%
+# Plotting multiple beachballs
+# ----------------------------
+#
+# Data of four earthquakes taken from USGS.
+# Provide lists.
+
+# Set up a dictionary
+aki_dict_multiple = {
+    "strike": [255, 173, 295, 318],
+    "dip": [70, 68, 79, 89],
+    "rake": [20, 83, -177, -179],
+    "magnitude": [7.0, 5.8, 6.0, 7.8],
+    "longitude": [-72.53, -79.61, 69.46, 37.01],
+    "latitude": [18.44, 0.90, 33.02, 37.23],
+    "depth": [13, 19, 4, 10],
+    "plot_longitude": [-70, -110, 100, 0],
+    "plot_latitude": [40, 10, 50, 55],
+    "event_name": [
+        "Haiti - 2010/01/12",
+        "Esmeraldas - 2022/03/27",
+        "Afghanistan - 2022/06/21",
+        "Syria/Turkey - 2023/02/06",
+    ],
+}
+# Convert to a pandas.DataFrame
+aki_df_multiple = pd.DataFrame(aki_dict_multiple)
+
+
+# %%
+# Adding a label
+# --------------
+#
+# Use the optional parameter ``event_name`` to add a label near the beachball,
+# e.g., event name or event date and time. Change the font of the the label text
+# by appending **+f** and the desired font (size,name,color) to the argument passed
+# to the ``scale`` parameter. Additionally, the location of the label relative to the
+# beachball [Default is ``"TC"``, i.e., Top Center]; can be changed by appending
+# **+j** and an offset can be applied by appending **+o** with values for *dx*\ /*dy*.
+# Add a colored [Default is white] box behind the label via the label ``labelbox``.
+# Force a fixed size of the beachball by appending **+m** to the argument passed to
+# the ``scale`` parameter.
+
+fig = pygmt.Figure()
+fig.coast(region="d", projection="N10c", land="lightgray", frame=True)
+
+fig.meca(spec=aki_df_multiple, scale="0.4c+m+f5p", labelbox="white@30", offset="+s0.1c")
+
+fig.show()
+
+
+# %%
+# Using size-coding and color-coding
+# ----------------------------------
+#
+# The beachball can be sized and colored by the quantities given as ``magnitude`` and
+# ``depth``, e.g., by moment magnitude or hypocentral depth, respectively. Use the
+# parameter ``cmap`` to pass the descired colormap. Now, the fills of the small circles
+# indicating the event locations are given by the colormap.
+
+fig = pygmt.Figure()
+fig.coast(region="d", projection="N10c", land="lightgray", frame=True)
+
+# Set up colormap and colorbar for hypocentral depth
+pygmt.makecpt(cmap="lajolla", series=[0, 20])
+fig.colorbar(frame=["x+lhypocentral depth", "y+lkm"])
+
+fig.meca(
+    spec=aki_df_multiple,
+    scale="0.4c+f5p",
+    offset="0.2p,gray30+s0.1c",
+    labelbox="white@30",
+    cmap=True,
+    outline="0.2p,gray30",
+)
+
+fig.show()
+
+# sphinx_gallery_thumbnail_number = 9