You can use GMT for this task but there are many more options in Julia.
The example below first grids the irregular points using the splines in tension method (Smith and Wessel, 1990). Once you have the grid, you can use your favorite plotting package to contour the data, including off course GMT.jl.
# 1 - INPUT DATA:
n = 200
xs, ys = 2π*(rand(n) .- 0.5), 2π*(rand(n) .- 0.5)
zs = 100*sin.(xs .* ys)
# 2 - GMT GRIDDING (Smith and Wessel [1990] Splines in tension):
using GMT
data = [xs ys zs]
x = y = LinRange(-π, π, 100)
G = GMT.surface(data, R=(extrema(x)..., extrema(y)...), inc=(step(x), step(y)), T=0.1)
