using GMT, GLVisualize

# ---------------------------------------------------------------------
function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T})
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

# If this doesn't work, dowload it manually
download("http://gmt.soest.hawaii.edu/projects/gmt/repository/show/trunk/test/potential/bathy_1m.nc")
fname = "bathy_1m.nc";
G = gmt("gmtread -Tg " * fname)

# Would be very nice if
# 1) Vectors of X and Y coordinates were accepted instead of obliging to provide meshgrids of them
# 2) Vectors of Float64 were accepted as well. When dealing with geographical coordinatres Float32 may be to crude.
xx,yy=meshgrid(convert(Array{Float32}, (G.x)), convert(Array{Float32}, (G.y)));

# Now, the result of this is kind of disastrous. My guess is that the visualization expects all three axes
# to be isometric but that's not the case here XX and XX are geographical coordinates an Z is in meters.
w,r = glscreen();
view(visualize(xx,yy,G.z, :surface))
r()

# The following example does not need any external file and works a bit (but just a bit) better.
# However, the origin is at the center of the figure instead of at lower left corner. Why? Because
# both xx and yy start at 0 instead of -w/2 (w = fig width)?

#=
G = gmt("surface -R0/150/0/100 -I1", rand(Float64,100,3)*150);
z = convert(Array{Float32}, (G.z-G.range[5])/(G.range[6]-G.range[5]));	# Normalize Z so we can see something.
xx,yy=meshgrid(convert(Array{Float32}, (G.x)), convert(Array{Float32}, (G.y)));
w,r = glscreen();
view(visualize(xx,yy,G.z, :surface))
r()
=#