Help me with cross section plot

Some body please help me using pygmt, projecting hypocenter to A-A1 line and to plot the cross-section along A-A1.
The following is my pygmt script and the resulting figure.
Thank…

import pygmt
import pandas as pd
data=pd.read_csv("G:/LAT_JUPYTER/PyGMT/sumatra_eq.csv")
fig = pygmt.Figure()
fig.basemap(region=[94, 108, -7, 7], projection="M8i", frame="a2g")
fig.coast(shorelines="1.5p,black")

fig.plot(
x=data.longitude,
y=data.latitude,
size=0.15*data.magnitude,
style="cc",
color="red",
pen="black")

fig.plot(x=[96,103],y=[-3,2],pen="2p,black")
fig.text(text="A", x=95.5, y=-3, font="24p,Helvetica-Bold,black")
fig.text(text="A1", x=103.2, y=2.5, font="24p,Helvetica-Bold,black")

pygmt.project(output="cross.dat",
x=data.longitude,
y=data.latitude,
z=data.depth,
center=["-C96/-3"],
endpoint=["-E103/2"],
width="-W-100/100")

fig.show()
fig.savefig("G:/LAT_JUPYTER/PyGMT/sumatra.pdf")
2 Likes

Hi @agusms! Welcome to the forum! So it will be a bit hard to know what’s happening without your data, but I think you want to try this for your pygmt.project function:

pygmt.project(
    output="cross.dat",
    x=data.longitude,
    y=data.latitude,
    z=data.depth,
    center="96/-3",
    endpoint="103/2",
    width="-100/100"
)

Specifically, you don’t need to add the -C/-E/-W, because those are what the center/endpoint/width aliases are for. After you have the “cross.dat” file, you will need to plot that projected data somehow using fig.plot. If you are converting from using a GMT script, you can post it here and we can help to translate if for you to PyGMT.

Oh, and if you need more help, it might be better to do a print(data) and provide the output, or upload your CSV as a ZIP file. That way we’ll know what things look like :smiley:

Hi @weiji14, thank you for your help. I attach zip file of GMT script, eq data zipzip.zip (188.4 KB) to plot sumatran seismicity and cross section of single line projection only.
Will you please help me translate it to PyGMT ? Is it possible also using pygmt.Figure.subplot() ?

Thanks

Hi @agusms, thanks for providing the script and data! The code is as below.

Original GMT cross.bat script:

rem To plot eq and cross-section

gmt makecpt -Cred,green,blue -T0,70,300,10000 > quakes.cpt

gmt begin sumatra pdf
  echo 100 -6 >line
  echo 104 -2 >>line
  echo 99.5 -6.2 A > word
  echo 104.5 -1.5 B >> word
  gmt coast -R95/115/-10/5 -JM4i -Bxa2g2 -Bya2g2pf3   -Glightbrown -B+t"Sumatra"  -W0.25p
  gmt plot eq.dat -Wfaint  -Sc0.15 -Cquakes.cpt
  gmt psxy  line   -JM  -W2
  gmt text word -JM -F+f15,Helvetica
  
  gmt project eq.dat -Q -C100/-6.0 -E104.0/-2 -Fpz -W-100/100 > cross.dat
  gmt basemap -JX10/-6 -R0/500/0/150 -Bxafg100+l"Distance"  -Byafg50+l"Depth" -BWSen -Y-8 
  gmt psxy cross.dat -JX  -Sc0.2 -W1 -Gred
  echo 10 10 A > word2
  echo 490 10 B >> word2
  gmt text word2 -JX -F+f13,Helvetica
  
gmt end

Translated PyGMT v0.5.0 script:

import pygmt


fig = pygmt.Figure()

pygmt.makecpt(cmap="red,green,blue", series="0,70,300,10000", output="quakes.cpt")

fig.coast(
    region=[95, 115, -10, 5],
    projection="M4i",
    frame=["xa2g2", "ya2g2pf3", '+t"Sumatra"'],
    land="lightbrown",
    shorelines="0.25p",
)
fig.plot(data="eq.dat", pen="faint", style="c0.15", cmap="quakes.cpt")
fig.plot(x=[100, 104], y=[-6, -2], projection="M", pen=2)
fig.text(x=99.5, y=-6.2, text="A", font="15,Helvetica")
fig.text(x=104.5, y=-1.5, text="B", font="15,Helvetica")

pygmt.project(
    data="eq.dat",
    unit=True,
    center=[100, -6.0],
    endpoint=[104.0, -2],
    convention="pz",
    width=[-100, 100],
    outfile="cross.dat",
)
fig.basemap(
    projection="X10/-6",
    region=[0, 500, 0, 150],
    frame=['xafg100+l"Distance"', 'yafg50+l"Depth"', "WSen"],
    yshift=-8,
)
fig.plot(data="cross.dat", projection="X", style="c0.2", pen=1, color="red")
fig.text(x=10, y=10, text="A", font="13,Helvetica")
fig.text(x=490, y=10, text="B", font="13,Helvetica")

fig.savefig(fname="sumatra.pdf")
fig.show()

produces:

Hopefully you find it helpful!

Yes, but since you are just using 2 panels, it might be easier just to use xshift/yshift or fig.shift_origin. If you’re interested in making more complicated subplots, have a look at the tutorial at https://www.pygmt.org/v0.5.0/tutorials/subplots.html.

Thank you very much @weiji14 for your help, it’s done.
Now I have to plot the contour of potential hazard at some region of Sumatra. If I am facing the diffulties, probably I need your help.

Thank you.
agus

1 Like

Hi, I have a bit similar problem. How to plot the cross section into north south and east west direction including the contour? something like the attached figure. Thankyou before
image

1 Like

I would approach this problem by using grdtrack to get the surface profile to plot in panels B and C.

… and build from bottom-to-top (at least for now):

Hi, do you have any example using grtrack?

if you do, could you share the code please?

Thank you so much

did you find any example to do this?

if you did, could you share how did you it?

Thank you so much.

Yes, there is an example in the docs for grdtrack - https://www.pygmt.org/latest/api/generated/pygmt.grdtrack.html#examples-using-pygmt-grdtrack.