Anybody please help me to make plot slab2 model with cross section on seismicity ?
I put my code and picture here
New folder.zip (660.9 KB)
import pygmt
fig = pygmt.Figure()
#Buat CPT
pygmt.makecpt(cmap="red,green,blue", series="0,70,300,10000", output="quakes.cpt")
fig.coast(
region=[115, 120, -10, -7],
projection="M4i",
frame=["xa1g0", "ya1g0pf3"],
land="226/226/226",
water='lightblue',
shorelines="0.25p",
)
fig.plot(data="eq.dat", pen="faint", style="c0.15", cmap="depth.cpt")
fig.plot(data="eq2.dat", pen=0.7, style="a0.6", cmap="depth.cpt")
fig.colorbar(
position="g115.2/-7.2+w3c/0.3c+h",
box=False,
# Set cmap to the "roma" CPT
cmap="depth.cpt",
# Label the x-axis "Velocity" and the y-axis "m/s"
frame=['xafg0+l" "'],
#frame=['x+250"Depth"'],
)
fig.plot(
region=[115, 120, -10, -7],
projection="M4i",
frame="a",
x=[117, 117, 118.5, 118.5, 117],
y=[-9.4, -7.5, -7.5, -9.4, -9.4],
pen="1p,black",
)
# cross section PNG triggered
#fig.plot(x=[117.852, 117.852], y=[-9.5, -7.5], projection="M", pen=1)
#fig.text(x=117.9, y=-9.6, text="A", font="15,Helvetica")
#fig.text(x=117.9, y=-7.4, text="B", font="15,Helvetica")
# cross section triggered
#fig.plot(x=[117.763, 117.763], y=[-9.5, -7.5], projection="M", pen=1)
#fig.text(x=117.65, y=-9.6, text="C", font="15,Helvetica")
#fig.text(x=117.65, y=-7.4, text="D", font="15,Helvetica")
pygmt.project(
data="eq.dat",
unit=True,
center=[117.852, -9.6],
endpoint=[117.852, -7.4],
convention="pz",
width=[-150, 150],
outfile="cross.dat",
)
pygmt.project(
data="eq2.dat",
unit=True,
center=[117.852, -9.6],
endpoint=[117.852, -7.4],
convention="pz",
width=[-100, 100],
outfile="cross2.dat",
)
fig.plot(
projection="X10/-6",
region=[0, 400, 0, 200],
frame=['xafg0+l"Distance"', 'yafg0+l"Depth"', "WSen"],
x=[0, 400],
y=[50, 50],
pen="0.5p,--",
yshift=-8,
)
fig.plot(x=[0, 400], y=[100, 100], pen="0.5p,--")
fig.plot(x=[0, 400], y=[150, 150], pen="0.5p,--")
fig.plot(data="cross.dat", projection="X", style="c0.2", pen=0.1, color="red")
fig.plot(data="cross2.dat", projection="X", style="a0.5", pen=0.5, color="166/221/238")
#fig.text(x=10, y=10, text="A", font="13,Helvetica")
#fig.text(x=490, y=10, text="B", font="13,Helvetica")
fig.show()
!