Slab2 model on cross section

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()

Hello Buha,

Look at the PyGMT example for plotting profiles of bathymetry.
https://www.pygmt.org/latest/gallery/images/cross_section.html#sphx-glr-gallery-images-cross-section-py

You can run a very similar profile through the Slab2 grid file.

++Eric