Swapped compression/extension fills in meca module with mt convention #3682
Description
Description of the problem
For certain combinations of strike (mt convention) of the beach ball by meca module gives wrongly swapped compression & expression fill. The following figure is a typical example of this problem, showing the mechanism solution with the rake angle changing by 1 degree.
The aki convention does not cause this problem.
The moment tensor in the above example has an equal plunges between P-and T-axis, so this issue of GMT might be relevant. But note that the problem also occurs in the mt convention, not in the principal axis representation (-Sy option in GMT or principal_axis convention in PyGMT).
I also noticed that we cound avoid the problem by adding a very small value (such as 0.0001 degree) to the rake angle.
Minimal Complete Verifiable Example
def sdr2moment(strike, dip, rake, M0):
"""calculate moment tensor component from strike, dip, and rake
Reference: Aki and Richards (2002) Eq (1) of Box 4.4
"""
sinδ = np.sin( np.deg2rad(dip ))
cosδ = np.cos( np.deg2rad(dip ))
sin2δ = np.sin(2 * np.deg2rad(dip ))
cos2δ = np.cos(2 * np.deg2rad(dip ))
sinλ = np.sin( np.deg2rad(rake ))
cosλ = np.cos( np.deg2rad(rake ))
sinφ = np.sin( np.deg2rad(strike))
cosφ = np.cos( np.deg2rad(strike))
sin2φ = np.sin(2 * np.deg2rad(strike))
cos2φ = np.cos(2 * np.deg2rad(strike))
Mxx = - M0 * (sinδ * cosλ * sin2φ + sin2δ * sinλ * sinφ * sinφ )
Mxy = M0 * (sinδ * cosλ * cos2φ + sin2δ * sinλ * sin2φ / 2 )
Mxz = - M0 * (cosδ * cosλ * cosφ + cos2δ * sinλ * sinφ )
Myy = M0 * (sinδ * cosλ * sin2φ - sin2δ * sinλ * cosφ * cosφ )
Myz = - M0 * (cosδ * cosλ * sinφ - cos2δ * sinλ * cosφ )
Mzz = M0 * ( sin2δ * sinλ )
return Mxx, Myy, Mzz, Myz, Mxz, Mxy
def meca_aki2mt(strike, dip, rake, Mw):
M0 = 10**( 1.5 * (Mw + 10.7) ) # in dyn-cm unit
exponent = int(np.floor(np.log10(M0)))
# strike, dip, rake -> moment tensor in the cartesian coordinate
Mxx, Myy, Mzz, Myz, Mxz, Mxy = sdr2moment(strike, dip, rake, M0/10**exponent)
# from cartesian to polar coordinates
Mrr, Mtt, Mff, Mtf, Mrf, Mrt = Mzz, Mxx, Myy, -Mxy, -Myz, Mxz
spec = {"mrr": Mrr, "mtt": Mtt, "mff": Mff,
"mrf": Mrf, "mrt": Mrt, "mtf": Mtf,
"exponent": exponent}
return spec
## main script
fig = pygmt.Figure()
strike = 75
dip = 30
rake = 180
fig.basemap(projection='X10c/10c', region=[0, 4, 0, 4], frame = 'wsen')
fig.text(x=2, y=3.8, font="12p,Helvetica,Red", text="aki")
fig.text(x=3, y=3.8, font="12p,Helvetica,Blue", text="mt")
fig.text(x=1, y=3.8, font="10p,Helvetica,Black", text="(@~f@~,@~d@~,@~l@~)", justify='RM')
for i in range(3):
rake = 179 + i
fig.text(x=1, y=3-i, text=f"({strike}, {dip}, {rake})", justify = 'RM')
fig.meca(spec={"strike": strike, "dip": dip, "rake": rake, 'magnitude': 5},
longitude= 2, latitude=3-i, scale="2c", compressionfill='red')
fig.meca(spec = meca_aki2mt(strike, dip, rake, 5),
longitude= 3, latitude=3-i, scale="2c", compressionfill='blue')
fig.plot(x=[2.5, 3.5, 3.5, 2.5, 2.5], y=[1.5, 1.5, 2.5, 2.5, 1.5], pen='thicker,orange')
fig.show()Full error message
No error message appeared.System information
PyGMT information:
version: v0.13.0
System information:
python: 3.12.5 | packaged by conda-forge | (main, Aug 8 2024, 18:36:51) [GCC 12.4.0]
executable: /home/tktmyd/miniforge3/envs/seismo24/bin/python
machine: Linux-6.8.0-45-generic-x86_64-with-glibc2.35
Dependency information:
numpy: 1.26.4
pandas: 2.2.2
xarray: 2024.9.0
netCDF4: 1.7.1
packaging: 24.1
contextily: None
geopandas: None
IPython: 8.27.0
rioxarray: None
gdal: 3.9.2
ghostscript: 10.03.1
GMT library information:
version: 6.5.0
padding: 2
share dir: /home/tktmyd/miniforge3/envs/seismo24/share/gmt
plugin dir: /home/tktmyd/miniforge3/envs/seismo24/lib/gmt/plugins
library path: /home/tktmyd/miniforge3/envs/seismo24/lib/libgmt.so
cores: 36
grid layout: rows
image layout:
binary version: 6.5.0