Stacking automatically generated cross-profiles
makecpt(cmap=:rainbow, range=(-5000,-2000))
grdimage("@spac_33.nc", shade=(azim=15, norm="e0.75"), proj=:merc, figsize=15)
# Select two points along the ridge
ridge_pts = [-111.6 -43.0; -113.3 -47.5];
# Plot ridge segment and end points
plot!(ridge_pts, region="@spac_33.nc", symbol=(symb=:circle, size=0.25),
fill=:blue, pen=(2,:blue))
# Generate cross-profiles 400 km long, spaced 10 km, samped every 2km
# and stack these using the median, write stacked profile
table = grdtrack("@spac_33.nc", ridge_pts, equidistant="400k/2k/10k+v", stack="m+sstack.txt")
plot!(table, pen=0.5)
# Show upper/lower values encountered as an envelope
env1 = gmtconvert("stack.txt", outcol="0,5")
env2 = gmtconvert("stack.txt", outcol="0,6", reverse=true, suppress=true)
env = [env1[1].data; env2[1].data]; # Concat the two matrices
plot!(env, region=(-200,200,-3500,-2000), proj=:linear, figsize=(15,7.5),
frame=(axes=:WSne, annot=:auto, ticks=:auto),
xaxis=(grid=1000, label="Distance from ridge (km)"),
ylabel="Depth (m)", fill=:lightgray, yshift=16)
plot!("stack.txt", pen=3)
text!(text_record([0 -2000], "MEDIAN STACKED PROFILE"), fill=:white, font=14,
justify=:TC, offset=(away=true, shift=0.25), show=true)
rm("stack.txt")
See also GMT ex33