HI,How can I use pygmt to extract the elevation curve between two points, as shown in the figure below

[Phd.Yin]

I want to plot a cross section and show the earthquakes.

Thanks!

[volkerk]

You can proceed as follows

1. Use gmt project to get coordinates along a specific track -> track.txt
2. Use track,txt in gmt grdtrack applied to a topography grid to get the respective elavation for the coordinate pairs which will add the interpolated elevation as a further collumn. -> track_elevation.txt
3. Plot third and fourth collumn of track_elevation.txt containing the distance on the track and elevation, respectively, e.g., using gmt plot.

[pwessel]

If it is just a straight profile from A to B you can bypass project and use grdtrack’s -E option to define the profile.

[Phd.Yin]

Thanks, but I ofen use the pygmt to plot. Is there some example by using pygmt?

[mdtanker]

Here’s a bit of code I use to do this, and the output figure. I hope it helps.

import pandas as pd
import numpy as np
import pygmt

# set start/stop points of line (here it's meters N/S in a polar stereographic projection), and grid to sample
a=(-590000,-1070000)
b=(-100000,-545000)
input_grid='tmp/BedMachine_surface_wgs_5k.nc'

# Create line from point a to point b with a row at every 1000m between lines.
df = pd.DataFrame(data = np.linspace(start = a, stop = b, num = 1000), columns = ["x", "y"])

# Make column of distance from point a
df['Distance'] = np.sqrt( (df.x-df.x.min())**2 + (df.y-df.y.min())**2 )

# Sample 'input_grid' at every 1000m along line from a to b, and save sampled data to new column
df = pygmt.grdtrack(points=df, grid=input_grid, newcolname=input_grid)

fig = pygmt.Figure()
fig.plot(x = df.Distance, y = df[input_grid], pen = '2p,red', frame = True)
fig.show()

[Phd.Yin]

Thanks ! That very helpful!