Hi - I’ve been using PyGMT to extract and plot seismic tomography models using some very simple algorithms. I’ve found code that work with PyGMT version 0.3.x no longer works with version 0.4.1
I get an error message from pygmt.grdinfo , where the grid metadata is not returned.
I suspect there is something very simple to do with the xarray dataset. Any help would be appreciated
HEre is my code with some lines commented to illustrate behaviors:
# read_plot_wang_et_al_2020_model.py
import xarray as xr
import pygmt
pygmt.show_versions()
# Wang et al (2020) seismic tomography model of Southern CA, available at
# http://ds.iris.edu/ds/products/emc-socalanat_vsrawang2020/
filename = "SoCal.ANAT-Vs+RA.Wang.2020.nc"
ds = xr.load_dataset(filename)
vsv = ds.vsv.isel(depth=25)
# Plot 2D section for the 25th depth level with xarray, works fine
vsv.plot()
# pygmt cannot get the grid file metadata
pygmt.grdinfo(vsv)
# Get the domain limits from the dataset
ds_domain = [
[ ds.longitude.values.min(), ds.longitude.values.min(), ds.longitude.values.max(), ds.longitude.values.max(), ds.longitude.values.min() ],
[ ds.latitude.values.min(), ds.latitude.values.max(), ds.latitude.values.max(), ds.latitude.values.min(), ds.latitude.values.min() ]
]
ds_domain
# Make map with PyGMT
fig = pygmt.Figure()
region_m = [ -122.0, -114.5, 32.0, 37.0 ]
fig.coast(
region=region_m,
projection="M6c",
land="grey",
water="lightblue",
shorelines=True,
frame="a",
)
pygmt.makecpt(cmap="seis", series="2.0/4.5/0.5", continuous=True, output="vs.cpt")
# PyGMT cannot get grid metadata
#fig.grdimage(grid=vsv, cmap="vs.cpt", nan_transparent=True)
#fig.grdimage(grid=vsv)
fig.plot(ds_domain[0], ds_domain[1])
fig.colorbar(cmap="vs.cpt", frame=["x+lVSV", "y+lkm/s"])
fig.show()
