Help me with cross section plot

Q&A PyGMT Q&A
1

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")
3 Likes
2

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:

3

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

4

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:

sumatra
sumatra1364×2115 201 KB

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 Making subplots — PyGMT.

2 Likes
5

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

2 Likes
6

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
7

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

8

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

9

Hi, do you have any example using grtrack?

if you do, could you share the code please?

Thank you so much

10

did you find any example to do this?

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

Thank you so much.

11

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

12

@agusms could you share the code, I attempted with the script provided on here but in not works for me. I also need to make a similar plot.

13

@Adnan I think this script could be useful for you. Here, pygmt.project is used to project the data.

14

@Esteban82 I attempted with this script, but it did not work well in my cases. I wonder what are the input for eq.dat file.

15

@Adnan I think that the data file is in this zip file.

16

Anyone who can help me with the above. The cross section is not cut the data.

import pygmt

fig = pygmt.Figure()

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

fig.coast(

    region=[25.0, 25.5, 35.0, 35.5],

    projection="M4i",

    frame=["xa2g2", "ya2g2pf3", '+t"Sumatra"'],

    land="lightbrown",

    shorelines="0.25p",

)

fig.plot(data="crosssection.dat", pen="faint", style="c0.15", cmap="quakes.cpt")

# Add a line connecting two points and label them as A and A'

fig.plot(x=[25.140133, 25.363331], y=[35.175001, 35.109202], projection="M", pen=2)

fig.text(x=25.140133, y=35.175001, text="A", font="15,Helvetica")

fig.text(x=25.363331, y=35.109202, text="A'", font="15,Helvetica")

pygmt.project(

    data="crosssection.dat",

    unit=True,

    center=[25.100133, 35.175001],

    endpoint=[25.363331, 35.109202],

    convention="pz",

    width=[-100, 100],

    outfile="cross_AA.dat",

)

fig.basemap(

    projection="X10/-6",

    region=[-100, 200, 0, 30],

    frame=['xafg100+l"Distance"', 'yafg50+l"Depth"', "WSen"],

    yshift=1,

)

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="arkaloxori.png")

fig.show()

crosssection.dat (30.3 KB)

17

can you plot the slab1 or slab2 in this cross section plot ?