Vertical Cross-Section of 3D tomographic data

Hello everyone,

I was able to create a horizontal slice of my 3D data (peakdelay_3_Degrees_Hz.txt) with the help from GMT community member. Now I want to create a vertical cross-section. I tried to re write the code to do so. My code produces a vertical slice along a latitude or longitude line (its far from perfect ) but I want to be able to choose the start point and end point and get a vertical slice along the line. I was thinking about pygmt.project but it does not work with 3D data if I am correct. May be some one could see my code and suggest me what I could change please ? I am also uploading output figure. Thanks!

import pygmt 
import pandas as pd
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
from mpl_toolkits.mplot3d import Axes3D

datafile = '/mnt/d/Uni/Master_Thesis/GMT/TXT_NLLOC/peakdelay_3_Degrees_Hz.txt'
df = pd.read_csv(datafile,delimiter=',', names=['longitude', 'latitude', 'depth', 'peakdelay', 'number'])#, 'a', 'b','c','d','e'])#,

lat_slice = 2.6

df_lat = df[df['latitude']==lat_slice]

lons0 = np.array(df_lat['longitude'])
lats0 = np.array(df_lat['latitude'])
depth0 = np.array(df_lat['depth'])
pd = np.array(df_lat['peakdelay'])

coordinates0 = np.column_stack((lons0, depth0))

minlon, maxlon = min(lons0), max(lons0)
mindep, maxdep = min(depth0), max(depth0)
step = 0.01

lons = np.arange(minlon, maxlon, step)
deps = np.arange(mindep, maxdep, step)

xintrp, yintrp = np.meshgrid(lons, deps)
z1 = griddata(coordinates0, pd, (xintrp, yintrp), method='cubic')#,  'cubic' 'linear' 'nearest'
xintrp = np.array(xintrp, dtype=np.float32)
yintrp = np.array(yintrp, dtype=np.float32)

da = xr.DataArray(z1, dims=('dep', 'long'), coords={'long': lons, 'dep': deps},)

frame = ['a1f0.25', 'WSen']

# color pallets
fig = pygmt.Figure()

lim=abs(max(-0.4, 0.4))#( #(pd.min(), pd.max())) | pd (-0.5, 0.5)

    cmap='blue,white,red', #'blue,white,red',

    region=[minlon, maxlon, mindep, maxdep],

fig.text(x=97.4, y=1, text="xxxx", font="6,Helvetica")

    frame='+l"log. Peak delay"' #

That is something grdinterpolate does but not sure if it has been wrapped in pyGMT yet.