The following works for w = w(x,y), but if you have only a numerical matrix w then I think you could define an interpolator function of (x,y) to be called:
using Contour, Plots; gr(dpi=600)
x = y = range(0, π, 100)
z = @. cos(x) + sin(y')
w(x,y) = @. cos(x)^2 - sin(y')^2
plot()
for cl in levels(contours(x,y,z))
for line in lines(cl)
xs, ys = coordinates(line) # coordinates of each contour line
plot!(xs, ys, st=:path, markershape=:none, line_z = w(xs,ys)) # use plotting package you prefer
end
end
Plots.current()
