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},)
pygmt.config(MAP_FRAME_TYPE="plain")
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)
pygmt.makecpt(
cmap='blue,white,red', #'blue,white,red',
series=f'-{lim}/{lim}/0.01',
continuous=True
)
fig.grdimage(
region=[minlon, maxlon, mindep, maxdep],
grid=da,
projection='M2i',
interpolation='l'
)
fig.text(x=97.4, y=1, text="xxxx", font="6,Helvetica")
fig.colorbar(
frame='+l"log. Peak delay"' #
)
fig.show()
