Contour plot for non-rectangular domain

[Gus_Hart]

There was a previous thread on this, but the use case was too different from mine to be informative.

I have a function defined over a triangle and would like to plot it as a contour plot but only in the triangular domain. Is there some way to do this? It’s not hard in mathematica (but seems not to translate to more “normal” plotting software).

For a concrete example to aim at, consider f(x,y) = cos(2pi x)*sin(2pi y) over the triangle with vertices (0,0), (1,0), (0,1). I’d like the plot to be blank outside the triangle.

[jacobadenbaum]

Plots.jl ignores NaN values, so while this may be a bit hacky, you can do something like this:

using Plots
f(x,y) = x < 1 - y ? cos(2 * pi * x) * sin(2 * pi * y) : NaN
grid   = LinRange(0,1, 1000)
contourf(grid, grid, f)

You can obviously generalize this to arbitrary domains so long as it is easy for you to check whether a point lies inside it.

[Edit: I misread the domain you were asking about in the initial response]

[Gus_Hart]

I like that. It works. If there was some way to define a plot region, that would be cleaner, but this totally does the trick. Thank you.

[rafael.guerra]

If the input data is a sparse grid instead of a function evaluated on a dense grid, the nice solution above using NaNs may leak:

Here below another approach using Dierckx.jl to interpolate a bounding polygon (and respective values) which is then merged with the input data and fed to PyPlot.jl’s tricontour().

The main advantages of this workflow is that irregular input data can be consumed and “arbitrary” polygons can be defined.

PS: Someone versed in PyPlot, may step-in with further advice.

using Dierckx, PyPlot

# Define some input data inside a triangle
f(x,y) = cos(2π*x) * sin(2π*y)
Np = 2_000; x = zeros(Np); y = zeros(Np)
for i in 1:Np
    x[i] = rand(0:0.01:1);  y[i] = rand(0:0.01:1-x[i])
end
z = f.(x,y)

# define bounding polygon vertices for contouring (here a triangle)
polygon = [0. 0.; 1. 0.; 0. 1.; 0. 0.]   # closed triangle (4 points)
N = size(polygon,1)

# Interpolate polygon border using linear parametric splines, and respective values using cubic
t = collect(1:N)
spl1 = ParametricSpline(t, polygon', k=1)  # linear splines
tfine = LinRange(0, N, 1000)  # polygon boundary samples = 1000
xypoly1 = evaluate(spl1,tfine)   # 2 x 1000 matrix
spl2 = Spline2D(x, y, z; kx=3, ky=3, s=1e-4)
zpoly2 = evaluate(spl2, xypoly1[1,:], xypoly1[2,:])  # cubic splines

# merge input data with polygon boundary
xnew = [x; xypoly1[1,:]]
ynew = [y; xypoly1[2,:]]
znew = [z; zpoly2]
xyz = [xnew ynew znew]
xyz = unique(xyz, dims= 1)   # keep only unique points

# plot the data support and filled contours
clf()
PyPlot.scatter(x,y,s=0.1)
plt.plot(xypoly1[1,:],xypoly1[2,:], c="red", lw=1)
clf()
PyPlot.tricontourf(xyz[:,1], xyz[:,2], xyz[:,3])
PyPlot.colorbar()
PyPlot.tricontour(xyz[:,1], xyz[:,2], xyz[:,3], colors="black", linewidths=0.5)

[joa-quim]

It was already possible to do this with GMT but it was a bit cumbersome, now with master version one can

f(x,y) = cos(2 * pi * x) * sin(2 * pi * y);
x = linspace(0,1,100);
G = mat2grid(f, x, x);           # because countourf cannot yet ingest functions

# Clip inside a triangle for simplicity but the shape could be anything and any number of shapes
contourf(G, clip=[0.1 0.1; 0.5 0.9; 0.9 0.1], fmt=:png, show=true)

[Gus_Hart]

What is GMT?

[joa-quim]

GMT.jl the best mapping package