Based on this stack overflow post, a non-linear invertible function can be used to define a secondary PyPlot axis scale. Two questions, if you will:
(1) Is it possible to do this with other plot backends?
(2) If instead of functions, we have arrays of monotonous data points providing the relation between the two axes, how to define interpolation functions that can be consumed by PyPlot as in example below?
NOTE: the simple example is just for illustration purposes
using PyPlot
# For medium with a velocity gradient: V(z) = k*z + V0
# Invertible function defines secondary axis in PyPlot plot
Depth2Time(depth) = @. log(k*depth + V0)/k - log(V0)/k
Time2Depth(time) = @. exp(k*(time + log(V0)/k))/k - V0/k
const k = 1.5; # 1/s
const V0 = 1500; # m/s
depth = 0:2:2000 # meter
yz = @. sin(2π*10*k*depth/V0) * exp(-100*k^2*(depth - 500)^2/V0^2)
fig, ax = plt.subplots()
fig.subplots_adjust(bottom=0.25)
ax.set_xlabel("Depth [m]", color=:blue)
ax.tick_params(axis="x", colors= "blue")
ax.set_ylabel("Amplitude")
ax.plot(depth, yz)
# THIS IS THE KEY PYPLOT LINE:
ax2 = ax.secondary_xaxis(-0.2, functions=(Depth2Time, Time2Depth), color=:red)
ax2.set_xlabel("Time [s]")
plt.show()
As explained in this matplotlib documentation, sometimes we want to relate the axes in a transform that is derived empirically. We can then set the forward and inverse transform functions to be linear interpolations between the two datasets.
However, PyPlot does not seem to consume Julia’s interpolation functions, issuing error:
RuntimeError: <PyCall.jlwrap (in a Julia function called from Python)
JULIA: MethodError: no method matching (::Interpolations.GriddedInterpolation{Float64, 1, Float64, Gridded{Linear}, Tuple{Vector{Float64}}})(::Matrix{Float64})
See MWE below:
using PyPlot, Interpolations
N = 200;
velocity0 = rand(1500:3000,N-1) # interval velocity [m/s]
depth0 = collect(LinRange(0,2,N)) # depth [km]
time0 = [0; cumsum(1000*diff(depth0)./velocity0)] # time [s]
itp_fwd = interpolate((time0,), depth0, Gridded(Linear()))
itp_inv = interpolate((depth0,), time0, Gridded(Linear()))
Time2Depth(time) = itp_fwd(time) # s to km
Depth2Time(depth) = itp_inv(depth) # km to s
time = collect(0:0.002: 0.8)
Depth2Time(Time2Depth(time)) ≈ time # checks that inverse transforms work fine
yt = @. sin(2π*30*time) * exp(-900*(time - 0.3)^2)
fig, ax = plt.subplots()
fig.subplots_adjust(bottom=0.25)
ax.set_xlabel("Time [s]", color=:blue)
ax.tick_params(axis="x", colors= "blue")
ax.set_ylabel("Amplitude")
ax.plot(time, yt)
plt.xlim(extrema(time))
# Follow key PyPlot command fails with functions from interpolations.jl:
ax2 = ax.secondary_xaxis(-0.2, functions=(Time2Depth, Depth2Time), color=:red)
ax2.set_xlabel("Depth [km]")
plt.show()
Any feedback would be welcome. Thanks.
The error message suggests that you just need to broadcast:
Time2Depth(time) = itp_fwd.(time) # s to km
Depth2Time(depth) = itp_inv.(depth) # km to s
Interpolators are implemented for scalars and (column) vectors, so one might get away with using vec(), but it seems safest to preserve the array layout.
1 Like
@Ralph_Smith, thank you very much indeed! 
I would not have figured this one out, as everything else was working fine without broadcasting those functions.
Brilliant.