Cross-sections with pygmt.project

Q&A PyGMT Q&A
1

Hi,

I am new to pygmt and would like to plot a cross-section with the new pygmt.project functionality. I can reproduce this example but instead of plotting earthquake locations, I would like to plot a 2-D section through a 3-D model, e.g. S362ANI using the S362ANI_kmps netCDF file, which can be downloaded here.

What would be the easiest way of passing the shear-wave velocities in S362ANI on to pygmt.project? I believe some interpolation is required as well?

pygmt.project(
    data="S362ANI_kmps.nc",
    unit=True,
    center=[100, -6.0], # A
    endpoint=[104.0, -2], # B
    convention="pz",
    width=[-100, 100],
    outfile="cross.nc",
)

will throw:
project [ERROR]: No data columns to work with - exiting

I tried to pass on an ASCII file with columns “longitude”, “latitude”, “depth”, “vs”, which does not throw an error but cross.nc is empty.

Regarding the plotting, I assume I will have to use pygmt.grdimage() using cross.nc instead of pygmt.plot() as done in the Sumatra earthquake location example linked above.

Really appreciate any help with this!

2

project does not handle grids. if you want a slice out of a 3-D netcdf cube then grdinterpolate is the one.