That definitely did the trick. Many thanks. I’m posting the code and the resulting plot here in case anyone else is interested in using it for their own plot. Hopefully it will save someone some time.
region =[-121.8,-121.6,45.27,45.43,-7000,3500]
perspective=[260,12]
fig = pygmt.Figure()
# Configuration for the 'current figure'.
pygmt.config(MAP_FRAME_TYPE="plain")
pygmt.config(FORMAT_GEO_MAP="ddd.xx")
pygmt.config(FORMAT_FLOAT_OUT='%.12lg')
cutgrd = pygmt.grdcut(grid=grd, region=region, projection="M8i")
# All seismicity
fig.plot3d(x=eqLonAll,
y=eqLatAll,
z=-1000*np.array(eqDepAll),
sizes=0.1 * (2**np.array(eqMagAll)),
style="uc",
color="gray",
zscale="0.001c",
perspective=perspective,
projection="M8i",
#frame=["xa", "yaf", "za", "wSnEZ"],
region=region)
# Summit swarm
fig.plot3d(x=eqLon,
y=eqLat,
z=-1000*np.array(eqDep),
sizes=0.1 * (2**np.array(eqMag)),
style="uc",
color="red",
zscale="0.001c",
perspective=perspective,
projection="M8i",
#frame=["xa", "yaf", "za", "wSnEZ"],
region=region)
# Flank swarm
fig.plot3d(x=eqLonF,
y=eqLatF,
z=-1000*np.array(eqDepF),
sizes=0.1 * (2**np.array(eqMagF)),
style="uc",
color="blue",
zscale="0.001c",
perspective=perspective,
projection="M8i",
#frame=["xa", "yaf", "za", "wSnEZ"],
region=region)
fig.grdview(
grid=cutgrd,
perspective=perspective,
#frame=["xa", "yaf", "za", "wSnE"],
projection="M8i",
zscale="0.001c",
# Set the surftype to "surface"
surftype="s",
shading="+a45",
# Set the CPT to "geo"
cmap="geo",
region=region,
transparency=20
)
fig.plot3d(x=stalon,
y=stalat,
z=staelev,
style="t0.5c",
color="black",
zscale="0.001c",
perspective=perspective,
projection="M8i",
#frame=["xa", "yaf", "za", "WsNeZ"],
frame=["za", "WsNeZ"],
region=region)
fig.show(method='external')