Feature request: contour that can handle curvilinear grids

[Gyoshi]

surface currently accepts curvilinear grids. It would be nice if contour and contourf does too for plotting of parameterised values.

contours would also use VertexSurfaceLike a la #748
Similar to #675.

For example, the pyplot backend via Plots.jl does this:

import Plots; Plots.pyplot()
r1, r2 = 1.5, 1/2
r(u,v) = ((r1 + r2*cos(v))*cos(u), (r1 + r2*cos(v))*sin(u), r1/2*sin(u)+r2*sin(v)) # torus
us = vs = range(0, 2pi, length=25)
xs, ys, zs = [[pt[i] for pt in r.(us, vs')] for i in 1:3]
Plots.contourf(xs, ys, zs)

[ffreyer]

Closed as duplicate of #742

[briochemc]

To track contours and heatmaps issues separately, what about reopening this issue while reducing the scope of #742 to heatmaps only? This issue is about contours only, while #742 is about contours and heatmaps, but almost everything discussed in #742 is about heatmaps.

[briochemc]

@ffreyer Sorry to bump it again, but I think this issue should be reopened: It remains unresolved and the scope is distinct from #742.

Surely I can't be the only one still wanting to plot contours for curvilinear grids with Makie, right? Unless there is a workaround that I don't know about?

[asinghvi17]

It seems like the standard approach is to compute the contour lines of the z matrix directly, and then transform the coordinates according to the curvilinear grid (using bilinear interpolation).

It seems easy enough to do. But it would be pretty slow.

see also:

• https://observablehq.com/@recifs/contours-on-a-curvilinear-grid
• Support curvilinear grids  d3/d3-contour#64

If you plot a contour or contourf into a GeoAxis or PolarAxis, this is exactly what happens already.

[briochemc]

Yes I get this impression but how bad is it? Is there an example in Julia that uses Interpolations.jl or another package in the wild? Is there a better solution? I'm not sure how the contour polygons are built but it seems to me that they could be built directly from the curvilinear x and y?

I also see some Delaunay / Voronoi triangulations being used... but this also goes over my head a bit.

Maybe also relevant:

• Makie's tricontourf
• TopoPlots.jl
• https://github.com/kakearney/contourfcmap-pkg/blob/master/contourfcmap/contourfcmap.m

[asinghvi17]

Yeah tricontourf is probably your best immediate bet for this. It may even give you continuous contour lines.

All of these approaches (with the exception of scattered interpolation based things like TopoPlots) seem to use this same bilinear interpolation idea. So it's very doable...just has to get done :D

[briochemc]

I think both tricontourf does the same thing as TopoPlots, no?

I still do feel like a package like Contour.jl that handles irregular x and y directly should be much faster and hopefully not so hard? My (uneducated) understadning from looking at this short article on marching squares (or Wikipedia) is that the contour/polygon points are positioned on a the center of mass on the grid lines (where mass is value of the data z).

[briochemc]

What about using Contour.jl instead of Isoband.jl (cc @jkrumbiegel)? Despite the reply to this issue, Contour.jl works with matrices of x and y and curvilinear grids! (as opposed to Isoband.jl) See this MWE:

using Contour, Test

y = -10:10
x = -20:20
# The curvilinear grid:
X = [x + 0.01y^3 for x in x, y in y]
Y = [y + 10cos(x/40) for x in x, y in y]

# Just a plot to illustrate:
xrows = vcat(([x; NaN] for x in eachrow(X))...)
xcols = vcat(([x; NaN] for x in eachcol(X))...)
yrows = vcat(([x; NaN] for x in eachrow(Y))...)
ycols = vcat(([x; NaN] for x in eachcol(Y))...)

# The data on that grid:
Z = cosd.(X) .* sind.(Y)
levs = -1:0.05:1
contourlevels = Contour.contours(X, Y, Z, levs)

using GLMakie
fig = Figure(size = (400, 400))

# curvilinear surface
ax = Axis(fig[1, 1])
srf = surface!(ax, X, Y, Z; shading = NoShading)
grd = lines!(ax, [xrows; xcols], [yrows; ycols], color=(:black, 0.2))
translate!(grd, 0, 0, 100) # draw the grid above the surface

# contour lines
for cl in levels(contourlevels)
    lvl = level(cl) # the z-value of this contour level
    for line in Contour.lines(cl)
        xs, ys = coordinates(line) # coordinates of this line segment
        ls = lines!(ax, xs, ys, color=:red) # pseuod-code; use whatever plotting package you prefer
        translate!(ls, 0, 0, 110 + lvl) # draw the contours above all
    end
end

fig

[briochemc]

BTW AFAIU Contour.jl does exactly what I mentioned just above and computes the polygons directly without any need for bilinear interpolation, so it should be fast, right?! But maybe there other issues with it? 🤷

[briochemc]

If I use the original example of this issue using poly instead of lines:

using Pkg
Pkg.activate(".")
using Contour, Test

r1, r2 = 1.5, 1/2
r(u,v) = ((r1 + r2*cos(v))*cos(u), (r1 + r2*cos(v))*sin(u), r1/2*sin(u)+r2*sin(v)) # torus
us = vs = range(0, 2pi, length=25)
X, Y, Z = [[pt[i] for pt in r.(us, vs')] for i in 1:3]

levs = -1.3:0.1:1.3
colorrange = extrema(levs)
colormap = :viridis
contourlevels = Contour.contours(X, Y, Z, levs)

using GLMakie
fig = Figure(size = (400, 400))

# curvilinear surface
ax = Axis(fig[1, 1])

# contour lines
for cl in levels(contourlevels)
    lvl = level(cl) # the z-value of this contour level
    for line in Contour.lines(cl)
        xs, ys = coordinates(line) # coordinates of this line segment
        ls = poly!(ax, xs, ys; color=lvl, colormap, colorrange)
        translate!(ls, 0, 0, 110 + lvl) # draw the contours above all
    end
end

cb = Colorbar(fig[2, 1], srf; vertical = false)

fig

I get

There are some issues with how to go from isolines to isobands, but this is an encouraging proof-of-concept, no?

[asinghvi17]

Contour.jl is also using the same approach :D

see https://github.com/JuliaGeometry/Contour.jl/blob/6d155559af262cd7e218d9a7a87b34bfe16bf9c1/src/interpolate.jl#L42-L63 - it traces using marching squares across cells, and then simply pops the transformed coordinates in.

It's trivial to enable matrix-like input for contour (I had no idea Contours.jl could do this, we should update the documentation!)

@eval Makie begin
# just overload the contourlines function for matrix like input
# nothing else needs to be done
function contourlines(x::AbstractMatrix{<: Real}, y::AbstractMatrix{<: Real}, z::AbstractMatrix{ET}, levels, level_colors, labels, T) where {ET}
    contours = Contours.contours(x, y, z, convert(Vector{ET}, levels))
    return contourlines(T, contours, level_colors, labels)
end

end

Using the torus example,

r1, r2 = 1.5, 1/2
r(u,v) = ((r1 + r2*cos(v))*cos(u), (r1 + r2*cos(v))*sin(u), r1/2*sin(u)+r2*sin(v)) # torus
us = vs = range(0, 2pi, length=25)
xs, ys, zs = [[pt[i] for pt in r.(us, vs')] for i in 1:3]
contour(xs, ys, zs)

This looks a bit wonky but it technically works. On your warped curvilinear grid example it's a bit better:

# issue from @briochemc - curvilinear grid

y = -10:10
x = -20:20
# The curvilinear grid:
X = [x + 0.01y^3 for x in x, y in y]
Y = [y + 10cos(x/40) for x in x, y in y]

# Just a plot to illustrate:
xrows = vcat(([x; NaN] for x in eachrow(X))...)
xcols = vcat(([x; NaN] for x in eachcol(X))...)
yrows = vcat(([x; NaN] for x in eachrow(Y))...)
ycols = vcat(([x; NaN] for x in eachcol(Y))...)

# The data on that grid:
# Z = cosd.(X) .* sind.(Y)
Z = hypot.(X, Y .- 10) # just to make it completely clear what this should look like
levs = -1:0.05:1
contourlevels = Makie.ContoursHygiene.Contour.contours(X, Y, Z, levs)

# plot
fig = Figure(size = (400, 400))
# curvilinear surface
ax = Axis(fig[1, 1])
srf = surface!(ax, X, Y, Z; shading = NoShading)
grd = lines!(ax, [xrows; xcols], [yrows; ycols], color=(:black, 0.2))
translate!(grd, 0, 0, 100) # draw the grid above the surface

contour!(ax, X, Y, Z; color = :red)

fig

[asinghvi17]

The patch for contourf is a bit more complicated - we have to go inside the function and create a completely new version.

For code cleanliness, we should probably abstract the contour polygon calculation to an external function, or patch Isoband.jl to do this directly. But for now I'm just copy-pasting the contourf plot function and changing some things around.

The main change here is that I am swapping out the x and y inputs to Isoband for the axes of z, and then using a linear interpolation to obtain the values of x and y.

function Makie.plot!(c::Contourf{<:Tuple{<:AbstractMatrix{<:Real}, <:AbstractMatrix{<:Real}, <:AbstractMatrix{<:Real}}})
    xs, ys, zs = c[1:3]

    c.attributes[:_computed_levels] = lift(c, zs, c.levels, c.mode) do zs, levels, mode
        Makie._get_isoband_levels(Val(mode), levels, vec(zs))
    end

    colorrange = lift(c, c._computed_levels) do levels
        minimum(levels), maximum(levels)
    end
    computed_colormap = lift(Makie.compute_contourf_colormap, c, c._computed_levels, c.colormap, c.extendlow,
                             c.extendhigh)
    c.attributes[:_computed_colormap] = computed_colormap

    lowcolor = Observable{RGBAf}()
    Makie.lift!(Makie.compute_lowcolor, c, lowcolor, c.extendlow, c.colormap)
    c.attributes[:_computed_extendlow] = lowcolor
    is_extended_low = lift(!isnothing, c, c.extendlow)

    highcolor = Observable{RGBAf}()
    Makie.lift!(Makie.compute_highcolor, c, highcolor, c.extendhigh, c.colormap)
    c.attributes[:_computed_extendhigh] = highcolor
    is_extended_high = lift(!isnothing, c, c.extendhigh)
    PolyType = typeof(Polygon(Point2f[], [Point2f[]]))

    polys = Observable(PolyType[])
    colors = Observable(Float64[])

    function calculate_polys(xs, ys, zs, levels::Vector{Float32}, is_extended_low, is_extended_high)
        empty!(polys[])
        empty!(colors[])

        levels = copy(levels)
        @assert issorted(levels)
        is_extended_low && pushfirst!(levels, -Inf)
        is_extended_high && push!(levels, Inf)
        lows = levels[1:end-1]
        highs = levels[2:end]

        # zs needs to be transposed to match rest of makie
        # NOTE: we are using the axes of `zs` here.
        isos = Makie.Isoband.isobands(axes(zs, 1), axes(zs, 2), zs', lows, highs)

        levelcenters = (highs .+ lows) ./ 2
        # NOTE: this is an interpolation of the coordinate grid
        point_interp = Makie.Interpolations.linear_interpolation(axes(xs), Point2f.(xs, ys))

        for (i, (center, group)) in enumerate(zip(levelcenters, isos))
            points = Point2f.(group.x, group.y)
            polygroups = Makie._group_polys(points, group.id)
            for rectilinear_polygroup in polygroups
                # NOTE: This is the only major change
                # we reproject the lines to curvilinear space
                polygroup = map(rectilinear_polygroup) do ring
                    map(ring) do point
                        point_interp(point[1], point[2])
                    end
                end
                outline = polygroup[1]
                holes = polygroup[2:end]
                push!(polys[], GeometryBasics.Polygon(outline, holes))
                # use contour level center value as color
                push!(colors[], center)
            end
        end
        polys[] = polys[]
        return
    end

    onany(calculate_polys, c, xs, ys, zs, c._computed_levels, is_extended_low, is_extended_high)
    # onany doesn't get called without a push, so we call
    # it on a first run!
    calculate_polys(xs[], ys[], zs[], c._computed_levels[], is_extended_low[], is_extended_high[])

    poly!(c,
        polys,
        colormap = c._computed_colormap,
        colorrange = colorrange,
        highclip = highcolor,
        lowclip = lowcolor,
        nan_color = c.nan_color,
        color = colors,
        strokewidth = 0,
        strokecolor = :transparent,
        shading = NoShading,
        inspectable = c.inspectable,
        transparency = c.transparency
    )
end

After this, we can try this in multiple scenarios:

1. Observable notebook from Fil

function warp_space(p::Point)
    return Point2(
        2 + cos(p[2]) - cos(p[1]),
        2 + sin(p[2]) + sin(p[1]),
        # hypot(p[2], p[1])
    )
end

cartesian_pointgrid = Point2f.(LinRange(0, 1, 100) .- 0.5, LinRange(0, 1, 50)' .- 0.5)
zs = splat(hypot).(cartesian_pointgrid) 
pointgrid = cartesian_pointgrid .|> warp_space

contourf(first.(pointgrid), last.(pointgrid), zs)

2. Torus example

r1, r2 = 1.5, 1/2
r(u,v) = ((r1 + r2*cos(v))*cos(u), (r1 + r2*cos(v))*sin(u), r1/2*sin(u)+r2*sin(v)) # torus
us = vs = range(0, 2pi, length=25)
xs, ys, zs = [[pt[i] for pt in r.(us, vs')] for i in 1:3]
contourf(xs, ys, zs)

3. Deformed grid example

# issue from @briochemc - curvilinear grid

y = -10:10
x = -20:20
# The curvilinear grid:
X = [x + 0.01y^3 for x in x, y in y]
Y = [y + 10cos(x/40) for x in x, y in y]

# Just a plot to illustrate:
xrows = vcat(([x; NaN] for x in eachrow(X))...)
xcols = vcat(([x; NaN] for x in eachcol(X))...)
yrows = vcat(([x; NaN] for x in eachrow(Y))...)
ycols = vcat(([x; NaN] for x in eachcol(Y))...)

# The data on that grid:
# Z = cosd.(X) .* sind.(Y)
Z = hypot.(X, Y .- 10) # just to make it completely clear what this should look like
levs = -1:0.05:1

# plot
fig = Figure(size = (400, 400))
# curvilinear surface
ax = Axis(fig[1, 1])
srf = surface!(ax, X, Y, Z; shading = NoShading)
grd = lines!(ax, [xrows; xcols], [yrows; ycols], color=(:black, 0.2))
translate!(grd, 0, 0, 100) # draw the grid above the surface

contourf!(ax, X, Y, Z)

fig

[asinghvi17]

If I don't end up getting to it feel free to PR these changes!