heatmap and contour plots with curvilinear x, y arguments?

[dingraha]

Matplotlib (and I think Matlab) allow a user to pass 2D arrays for the x and y coordinates when creating a contour or filled contour plot. This is handy when the plot domain isn't rectangular. Makie, however, requires the x and y coordinates to be 1D arrays, which restricts the plot to a rectangular domain. Would it be possible to have the same 2D array functionality in Makie?

Thanks!

Daniel

[SimonDanisch]

See: https://github.com/JuliaPlots/Makie.jl/issues/675, https://github.com/JuliaPlots/Makie.jl/issues/720.

[dingraha]

Excellent, thanks. I hadn't found those open issues.

[briochemc]

I would recommend reopening this issue because I think it is a step beyond #675. Taking this screenshot from Wikipedia to illustrate what I mean (and sorry it should be 2D images instead of 3D),

My understanding is that heatmap currently:

• handles the cartesian grids (left) because it only takes in the extrema of the supplied x and y
• does not handle irregularly spaced rectilinear grids (middle), the feature requested in #675
• does not handle curvilinear grids (right), which is what this issue is about if I understand correctly.

[jdeldre]

I'd like to bump this, and stress particularly the desire for contour plots on 2-d curvilinear grids, as shown in the right image in @briochemc's post above. The underlying x,y points are supplied as arrays of the same size as the field data to be plotted, but a key point is that they can be smoothly mapped from a Cartesian grid. I have used PyPlot for years to access the matplotlib contour function mentioned by @dingraha, but would prefer to do this is Makie.jl.

[briochemc]

Gentle bump: Can we reopen this issue?

[asinghvi17]

Just to clarify, the missing feature here is support for arbitrary grids (a superset of rectilinear and curvilinear), right? Rectilinear grids are already supported as I understand.

[briochemc]

I reckon curvilinear is the superset but yes!

(The superset of curvilinear would be unstructured meshes but these are not indexed by a 2D array.)

[asinghvi17]

So it turns out that if you just construct a mesh in which no faces share any vertices, and assign colors appropriately, you can get a face-colored mesh. I figured out a way to convert from a mesh and a vector of face colors to a mesh with effectively per-face colors, via https://stackoverflow.com/questions/1664402/how-can-i-specify-per-face-colors-when-using-indexed-vertex-arrays-in-opengl-3-x.

The code for that was basically:

rs = LinRange(0, 1, 10)
θs = LinRange(0, π/2, 10)

polar2cart(r, θ) = Point2f(r * sin(θ), r * cos(θ))

pointgrid = polar2cart.(rs, θs')

vertices, faces, colors = curvilinear_grid_mesh(first.(pointgrid), last.(pointgrid), zeros(size(pointgrid)), rand(RGBAf, 10, 10))

mesh(vertices, faces; color = colors, shading = false, axis = (aspect = DataAspect(),))

so until someone gets per-face coloring into GLMakie as a possibility, you can use the code below.

Code for curvilinear_grid_mesh

```julia """ curvilinear_grid_mesh(xs, ys, zs, colors) Tesselates the grid defined by `xs` and `ys` in order to form a mesh with per-face coloring given by `colors`. The grid defined by `xs` and `ys` must have dimensions `(nx, ny) == size(colors) .+ 1`, as is the case for heatmap/image. """ function curvilinear_grid_mesh(xs, ys, zs, colors = zs) nx, ny = size(zs) ni, nj = size(colors) @assert (nx == ni+1) & (ny == nj+1) "Expected nx, ny = ni+1, nj+1; got nx=$nx, ny=$ny, ni=$ni, nj=$nj. nx/y are size(zs), ni/j are size(colors)." input_points_vec = Makie.matrix_grid(identity, xs, ys, zs) input_points = reshape(input_points_vec, size(colors) .+ 1) triangle_faces = Vector{GeometryBasics.TriangleFace{UInt32}}() triangle_points = Vector{Point3f}() triangle_colors = Vector{eltype(colors)}() sizehint!(triangle_faces, size(input_points, 1) * size(input_points, 2) * 2) sizehint!(triangle_points, size(input_points, 1) * size(input_points, 2) * 2 * 3) sizehint!(triangle_colors, size(input_points, 1) * size(input_points, 2) * 3) point_ind = 1 @inbounds for i in 1:(size(colors, 1)) for j in 1:(size(colors, 2)) # push two triangles to make a square # first triangle push!(triangle_points, input_points[i, j]) push!(triangle_points, input_points[i+1, j]) push!(triangle_points, input_points[i+1, j+1]) push!(triangle_colors, colors[i, j]); push!(triangle_colors, colors[i, j]); push!(triangle_colors, colors[i, j]) push!(triangle_faces, GeometryBasics.TriangleFace{UInt32}((point_ind, point_ind+1, point_ind+2))) point_ind += 3 # second triangle push!(triangle_points, input_points[i+1, j+1]) push!(triangle_points, input_points[i, j+1]) push!(triangle_points, input_points[i, j]) push!(triangle_colors, colors[i, j]); push!(triangle_colors, colors[i, j]); push!(triangle_colors, colors[i, j]) push!(triangle_faces, GeometryBasics.TriangleFace{UInt32}((point_ind, point_ind+1, point_ind+2))) point_ind += 3 end end return triangle_points, triangle_faces, triangle_colors end ```

Would people like to see this as a recipe?

Here is some code which takes a mesh and a vector of colors, and returns points, faces and colors to plot a per-face colored mesh.

Code for color_per_face

```julia function color_per_face(msh::GeometryBasics.Mesh{N, T, FaceType}, color::AbstractVector{ColorType}) where {N, T, FaceType, ColorType} positions = msh.position faces = GeometryBasics.faces(msh) @show N T FaceType facevec = Vector{GeometryBasics.TriangleFace{Int64}}() # facemat = zeros(UInt32, length(faces), 3) posvec = Vector{GeometryBasics.Point{N, T}}() colvec = Vector{ColorType}() sizehint!(facevec, length(faces)) sizehint!(posvec, length(faces) * 3) sizehint!(colvec, length(faces) * 3) pos_ind = 1 face_ind = 1 @inbounds for i in 1:length(faces) v1, v2, v3 = GeometryBasics.value.(faces[i].data) push!(posvec, positions[v1]) push!(posvec, positions[v2]) push!(posvec, positions[v3]) push!(colvec, color[i]) push!(colvec, color[i]) push!(colvec, color[i]) push!(facevec, GeometryBasics.TriangleFace((pos_ind, pos_ind+1, pos_ind+2))) pos_ind += 3 # 3 elements added to the vector so push 3 slides. end return posvec, facevec, colvec end newpos, newface, newcolors = color_per_face(msh, 1:length(faces(msh))) mesh(newpos, newface; color = newcolors, shading = false) ```

[SimonDanisch]

Yeah I don't think there will be per face coloring (since you can't do it in OpenGL) , so this is the way;)

[asinghvi17]

I guess that the most easily interpretable way would be a recipe called facecolormesh? Then, we could have an IrregularHeatmap which basically acts like a discretized version of surface (including 3d support if desired).

[jkrumbiegel]

Isn't it weird though that you effectively need four times the vertices because each face needs a duplication? In matlab it seems that per face coloring is the default, and their wireframes also look better with that because they are part of the primitive. Wonder how close we can get to that.

[SimonDanisch]

Not sure what matlab does, but relatively sure that they'll need to duplicate vertices for that. Wireframes are surprisingly hard, and the best quality can only be achieved with a special purpose shader written only for that use case.

[SimonDanisch]

Although I guess one idea would be to use a TextureBuffer to look up the color via the face index in the shader... Would be interesting to try that out.

[jkrumbiegel]

That's how I thought it would work, but no idea about opengl internals

[briochemc]

FWIW, data on curvilinear earth grids can come with the vertices of each face as metadata (e.g., as two (4,m,n)-sized arrays, one for longitude of vertices and one for latitudes). For example, this is the 3D array of longitude vertices for one of the CMIP6 models with a curvilinear grid:

julia> areads["lon_vertices"]
lon_vertices (4 × 360 × 300)
  Datatype:    Float32
  Dimensions:  vertices × i × j
  Attributes:
   units                = degrees_east

julia> areads["lon_vertices"][:,:,:]
4×360×300 Array{Float32, 3}:
[:, :, 1] =
 80.0  81.0  82.0  …  78.0  79.0
 81.0  82.0  83.0     79.0  80.0
 81.0  82.0  83.0     79.0  80.0
 80.0  81.0  82.0     78.0  79.0

[:, :, 2] =
 80.0  81.0  82.0  …  78.0  79.0
 81.0  82.0  83.0     79.0  80.0
 81.0  82.0  83.0     79.0  80.0
 80.0  81.0  82.0     78.0  79.0

;;; …

[:, :, 299] =
 80.0     80.0709  …  79.9291
 80.0709  80.1418     80.0
 80.0355  80.071      80.0
 80.0     80.0355     79.9645

[:, :, 300] =
 80.0     80.0355  …  79.9645
 80.0355  80.071      80.0
 80.0     80.0        80.0
 80.0     80.0        80.0

This type of curvilinear "tripolar" grid is fairly common for ocean models (looks like this image below) to avoid singularities at the North pole / Arctic ocean):

[asinghvi17]

Huh, interesting! I'd be curious to see whether this sort of grid would work with the function I posted - could you try it out, or post some data which I can use for a test in GeoMakie? Otherwise, I posted another function which you could probably adapt to this usecase.

[briochemc]

I'm not sure how to try it (only had time to skim your code) but this MWE will download the data you can try for yourself:

using Zarr, YAXArrays, Dates
using CSV
using DataFrames, DataFramesMeta

# Table of runs with google storage URLs
CMIP6_stores = CSV.read(download("https://storage.googleapis.com/cmip6/cmip6-zarr-consolidated-stores.csv"), DataFrame)

# Select a variable/model run with a curvilinear grid
stores = @chain CMIP6_stores begin
    @subset @byrow :variable_id=="areacello"
    @subset @byrow :table_id=="Ofx"
    @subset @byrow :source_id=="ACCESS-CM2"
    @subset @byrow :experiment_id=="historical"ensembles?)
    @subset @byrow :member_id=="r1i1p1f1"
end
store = only(stores.zstore)

# Open with Zarr
g = open_dataset(zopen(store, consolidated=true))

# Data to plot
lonv = g["vertices_longitude"][:,:,:]
latv = g["vertices_latitude"][:,:,:]
area2D = g["areacello"][:,:]

It should give you the following arrays:

lonv = g["vertices_longitude"][:,:,:]
4×360×300 reshape(::Array{Union{Missing, Float64}, 3}, 4, 360, 300) with eltype Union{Missing, Float64}:
[:, :, 1] =
 80.0  81.0  …  78.0  79.0
 81.0  82.0     79.0  80.0
 81.0  82.0     79.0  80.0
 80.0  81.0     78.0  79.0

[:, :, 2] =
 80.0  81.0  …  78.0  79.0
 81.0  82.0     79.0  80.0
 81.0  82.0     79.0  80.0
 80.0  81.0     78.0  79.0

;;; …

[:, :, 299] =
 80.0     …  79.9291
 80.0709     80.0
 80.0355     80.0
 80.0        79.9645

[:, :, 300] =
 80.0     …  79.9645
 80.0355     80.0
 80.0        80.0
 80.0        80.0

julia> latv = g["vertices_latitude"][:,:,:]
4×360×300 reshape(::Array{Union{Missing, Float64}, 3}, 4, 360, 300) with eltype Union{Missing, Float64}:
[:, :, 1] =
 -78.0     …  -78.0
 -78.0        -78.0
 -77.7532     -77.7532
 -77.7532     -77.7532

[:, :, 2] =
 -77.7532  …  -77.7532
 -77.7532     -77.7532
 -77.5057     -77.5057
 -77.5057     -77.5057

;;; …

[:, :, 299] =
 65.0     …  65.4183
 65.4183     65.0
 65.4191     65.0
 65.0        65.4191

[:, :, 300] =
 65.0     …  65.4191
 65.4191     65.0
 65.4193     65.0
 65.0        65.4193

julia> area2D = g["areacello"][:,:]
360×300 Matrix{Float32}:
 0.0  0.0  0.0  …  0.0  0.0
 0.0  0.0  0.0     0.0  0.0
 ⋮              ⋱
 0.0  0.0  0.0     0.0  0.0
 0.0  0.0  0.0     0.0  0.0

[asinghvi17]

Got the code to work by removing the ensembles? in the DataFramesMeta block.

I'm still not 100% sure how the data is structured. Is it a quad mesh with each quad defined by Point2.(lonv[:, i, j], latv[:, i, j])?

lines(vec(lonv), vec(latv); linewidth=0.1) heatmap(area2D)

[briochemc]

Ah yes sorry about the typo (failed to properly erase a comment)...

Yes you are correct! You might have to be a bit careful around the 0°/360° cutoff so that quads don't jump from one side to the other but yes these plots look correct to me.

[briochemc]

If I skip the cells with zero area, I get this sort of thing (again the quads that jump over the 0°/360° cutoff are an issue):

[asinghvi17]

Couple of problems but overall seems to work...

using Zarr, YAXArrays, Dates
using CSV
using DataFrames, DataFramesMeta

# Table of runs with google storage URLs
CMIP6_stores = CSV.read(Downloads.download("https://storage.googleapis.com/cmip6/cmip6-zarr-consolidated-stores.csv"), DataFrame)

# Select a variable/model run with a curvilinear grid
stores = @chain CMIP6_stores begin
    @subset @byrow :variable_id=="areacello"
    @subset @byrow :table_id=="Ofx"
    @subset @byrow :source_id=="ACCESS-CM2"
    @subset @byrow :experiment_id=="historical"
    @subset @byrow :member_id=="r1i1p1f1"
end
store = only(stores.zstore)

# Open with Zarr
g = open_dataset(zopen(store, consolidated=true))

# Data to plot
lonv = g["vertices_longitude"][:,:,1:300] .|> Float64
latv = g["vertices_latitude"][:,:,1:300] .|> Float64
area2D = g["areacello"][:,1:300] .|> Float64

quad_points = vcat([Point2{Float64}.(lonv[:, i, j], latv[:, i, j]) for i in 1:size(lonv, 2), j in 1:size(lonv, 3)]...)
quad_faces = vcat([begin; j = (i-1) * 4 + 1; [j j+1  j+2; j+2 j+3 j]; end for i in 1:length(quad_points)÷4]...)
colors_per_point = vcat(fill.(vec(area2D), 4)...)

GLMakie.activate!()
mesh(quad_points, quad_faces; color = colors_per_point, shading = false)

[briochemc]

Ha great, that's neat! Here's the same with some fixes for longitudes

(For some reason I cannot upload any images anymore on GitHub, but the simple mesh plot without GeoMakie is beautiful I promise! The code below's first fig should produce it for you.)

Although there is still something wrong near the north pole with GeoMakie for quads that go over the 180° line:

(using max(min(lonv, 180), -180) suffices to "fix" these north pole issues — not shown because cannot upload anymore...)

using Zarr, YAXArrays, Dates
using GLMakie, GeoMakie
using CairoMakie.GeometryBasics
using CSV
using DataFrames
using DataFramesMeta

@time CMIP6_stores = CSV.read(download("https://storage.googleapis.com/cmip6/cmip6-zarr-consolidated-stores.csv"), DataFrame)

stores = @chain CMIP6_stores begin
    @subset @byrow :variable_id=="areacello"
    @subset @byrow :table_id=="Ofx"
    @subset @byrow :source_id=="ACCESS-CM2"
    @subset @byrow :experiment_id=="historical" 
    @subset @byrow :member_id=="r1i1p1f1"
end

store = only(stores.zstore)
g = open_dataset(zopen(store, consolidated=true))
lonv = g["vertices_longitude"][:,:,:]
# make sure quads are not too distorted
loninsamewindow(l1, l2) = mod(l1 - l2 + 180, 360) + l2 - 180
lon = mod.(g["longitude"][:,:] .+ 180, 360) .- 180
lonv = loninsamewindow.(lonv, reshape(lon, (1, size(lon)...)))
latv = g["vertices_latitude"][:,:,:]
area2D = replace(g["areacello"][:,:], 0.0=>NaN)

quad_points = vcat([Point2{Float64}.(lonv[:, i, j], latv[:, i, j]) for i in 1:size(lonv, 2), j in 1:size(lonv, 3)]...)
quad_faces = vcat([begin; j = (i-1) * 4 + 1; [j j+1  j+2; j+2 j+3 j]; end for i in 1:length(quad_points)÷4]...)
colors_per_point = vcat(fill.(vec(area2D), 4)...)

GLMakie.activate!()
fig = Figure()
ax = Axis(fig[1,1])
mesh!(ax, quad_points, quad_faces; color = colors_per_point, shading = false)
xlims!(ax, (-180, 180))
ylims!(ax, (-90, 90))

fig = Figure(resolution = (1200,600))
ax = GeoAxis(fig[1,1]; dest = "+proj=moll")
mesh!(ax, quad_points, quad_faces; color = colors_per_point, shading = false)
cl=lines!(ax, GeoMakie.coastlines(), color = :black, linewidth=0.85)
translate!(cl, 0, 0, 1000)
fig

[asinghvi17]

Neat! I guess that's the last stretch to go on for this - getting geo grids adapted for GeoMakie.

(Sidenote: a slightly shorter way is clamp.(lonv. -180, 180))

Do you mind if I use this data for an example in GeoMakie?

[briochemc]

Yes, please go ahead! (It's not my data it's CMIP's, but it's public and great if there are examples on how to use it out there!)

[briochemc]

BTW here is the plot I could not upload last time around, for completeness.

My fix was to shift the longitudes of the quad vertices. I.e., if the quad center had longitude l2, then I applied

# make sure quads are not too distorted
loninsamewindow(l1, l2) = mod(l1 - l2 + 180, 360) + l2 - 180

to ensure that the vertex longitude l1 is within (l2 - 180, l2 + 180). But you might want to address the wrapping of cells around the edges of the projection differently because I'm guessing my bandaid fix will fail for some projections. FWIW, I think this PR is the pcolormesh fix in Cartopy (not sure what they do exactly).

[asinghvi17]

Yes it would probably need to be more general for e.g. interrupted Good homolosine and friends, as well as that Clima square projection...

Here's the tripolar mesh in 3D! The method is kinda hacky though, we need better ways of figuring out limits for these kinds of situations.

(and script, adding on to yours):

using Geodesy
transf = ECEFfromLLA(WGS84())

mesh_transform = Makie.PointTrans{3}() do p::Point3
    return Point3f((transf(LLA(p[2], p[1], p[3])) ./ 5f4)...) # use a scale factor to avoid Float32 inaccuracy
end 

f, a, p = mesh(quad_points, quad_faces; color = colors_per_point, axis = (type = LScene,))
p.transformation.transform_func[] = mesh_transform
p.shading = false
f

cc @rafaqz - you could probably do something similar, possibly with surface, if you can decompose Tyler plots to a recipe instead of the struct form it has now. The same would go for vector tile map plots.

[briochemc]

Just noticed this related post on discourse which made me think to ping @juliohm here since this seems to fall within the scope of Meshes.jl recipes AFAIU. Not sure how Meshes.jl work for spherical data projections.

[juliohm]

The recipes in Meshes.jl can plot all sorts of discretizations over the sphere, assigning colors to either elements or vertices of the mesh as explained in the linked post. Here is an example of a discretization with m x n meridians and parallels:

using Meshes
import GLMakie as Mke

s = Sphere((0,0,0))
m = discretize(s, RegularDiscretization(50, 30))

Mke.plot(m, color = 1:nelements(m), showfacets = true, facetcolor = :teal)

Projecting the geographic coordinates to projected coordinates is a simple task. We just need an additional design change in Meshes.jl to plug these projections seamlessly in the large list of transforms that are already available. This will be better explained in the book.

[briochemc]

@juliohm thanks!

Maybe I am missing something but ngons that cross an edge of the projection do not seem trivial to deal with (to me). For example, see how these curvilinear-grid quadrangles that cross the 180° meridian mess up the North pole.

My ugly band-aid solution in the zoomed in image below was to simply mask out those edge-crossing culprits:

I think ideally the crossing ngons should be split into smaller ngons that do not cross the edge for them to be properly displayed in the projection, but surely that's not that trivial, right?

[juliohm]

Agree it is not that simple without the machinery in Meshes.jl. There we have access to connectivity information, can disconnect, remove, insert quadrangles in the mesh as necessary to accommodate topological issues.

Even though I'm not planning to cover these low-level details in the book (targeting a different audience), the contents will certainly clarify the vision we have for the future.

I will share the draft with you as soon as it is ready. Your feedback is appreciated.

Em qui., 21 de set. de 2023 22:13, Benoit Pasquier @.***> escreveu:

@juliohm https://github.com/juliohm thanks!

Maybe I am missing something but ngons that cross an edge of the projection do not seem trivial to deal with (to me). For example, see how these curvilinear-grid quadrangles that cross the 180° meridian mess up the North pole. [image: Screenshot 2023-09-22 at 11 08 15 am] https://user-images.githubusercontent.com/4486578/269801282-d22b44e1-a8a0-46e3-9334-6ba462701179.png

My ugly band-aid solution in the zoomed in image below was to simply mask out those edge-crossing culprits: [image: Screenshot 2023-09-22 at 11 06 01 am] https://user-images.githubusercontent.com/4486578/269800950-90d162be-49c7-4846-8e0a-756c317ecc75.png

I think ideally the crossing ngons should be split into smaller ngons that do not cross the edge for them to be properly displayed in the projection, but surely that's not that trivial, right?

— Reply to this email directly, view it on GitHub https://github.com/MakieOrg/Makie.jl/issues/742#issuecomment-1730593969, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZQW3LR6SPPTFJKNOUAXBLX3TQ2BANCNFSM4TA3FQJA . You are receiving this because you were mentioned.Message ID: @.***>

[briochemc]

I'll ask a probably dumb question then: Should GeoMakie.jl depend on Meshes.jl to leverage these tools?

[juliohm]

I don't think it is a dumb question @briochemc. You are asking what many people have asked in the past.

My understanding is that GeoMakie.jl creates the appropriate axis for visualization of projected coordinates, so that whenever we send the projected coordinates to it, it will display in the right place on the screen. Similar to what the Makie.PolarAxis does when we pass coordinates that are polar: https://juliaearth.github.io/GeoStatsDocs/stable/variography/empirical.html#Varioplanes

The process of converting coordinates and manipulating the geospatial domain, however, is something that can live anywhere, and we are building a very sophisticated stack for that. So in theory you will be able to create your GeoMakie axis and call our recipe on it to get the correct visualization, our recipe will take care of handling the topological issues.

[briochemc]

Gentle bump: Any updates on this issue?

[SimonDanisch]

On which of the topics being discussed here?

[briochemc]

🙏 Sorry, I should have been more specific. I was talking about the original issue: heatmap and contour plots with curvilinear x, y arguments (passed as matrices, i.e., 2D arrays). It looks like all the ingredients for heatmap have been sprinkled around here but I'm not sure if someone put it all together somewhere else and could submit a PR to finally close this, and I am unaware of any attempt with contour, hence my gentle-bump question 😃

Also, apologies for polluting the thread with slightly-off-topic discussions about the specific issues related to earth projections that probably would be better suited to discussions in GeoMakie.

[SimonDanisch]

does surface(xmat2d, ymat2d, zeros(size(zmat2d)), color=zmat2d, shading=NoShading, interpolate=false) work for that?

[briochemc]

does surface(xmat2d, ymat2d, zeros(size(zmat2d)), color=zmat2d, shading=NoShading, interpolate=false) work for that?

That throws for me (no interpolate kwarg):

The available plot attributes for Surface{Tuple{Matrix{Float64}, Matrix{Float64}, Matrix{Float32}}} are:

alpha      color     colorrange  depth_shift  fxaa      inspectable      inspector_hover  invert_normals  material  nan_color  shading    space     ssao            transparency
backlight  colormap  colorscale  diffuse      highclip  inspector_clear  inspector_label  lowclip         model     overdraw   shininess  specular  transformation  visible

Generic attributes are:

cycle  dim_conversions  label  model  rasterize  transformation  xautolimits  yautolimits  zautolimits

Stacktrace:
 [1] validate_attribute_keys(plot::Surface{Tuple{Matrix{Float64}, Matrix{Float64}, Matrix{Float32}}})
   @ MakieCore ~/.julia/packages/MakieCore/f3yyf/src/recipes.jl:733
 [2] push!(scene::Scene, plot::Plot)
   @ Makie ~/.julia/packages/Makie/rEu75/src/scenes.jl:459
 [3] plot!
   @ ~/.julia/packages/Makie/rEu75/src/interfaces.jl:390 [inlined]
 [4] plot!(ax::Axis, plot::Surface{Tuple{Matrix{Float64}, Matrix{Float64}, Matrix{Float32}}})
   @ Makie ~/.julia/packages/Makie/rEu75/src/figureplotting.jl:412
 [5] _create_plot!(::Function, ::Dict{Symbol, Any}, ::Axis, ::Matrix{Float64}, ::Vararg{Matrix{Float64}})
   @ Makie ~/.julia/packages/Makie/rEu75/src/figureplotting.jl:381
 [6] surface!(::Axis, ::Vararg{Any}; kw::@Kwargs{color::Matrix{Float64}, shading::MakieCore.ShadingAlgorithm, interpolate::Bool})
   @ MakieCore ~/.julia/packages/MakieCore/f3yyf/src/recipes.jl:440
 [7] top-level scope
   @ REPL[25]:1

[SimonDanisch]

Ah, that change is still not tagged - it was released quite a while ago, but the attribute erroring refactor overlooked it 😅 You try master ;)

[briochemc]

Just a heads up. This now works in GLMakie v0.10.9:

julia> x = 1:3
1:3

julia> y = 1:4
1:4

julia> X = repeat(x, outer=(1, length(y)))
^[[B3×4 Matrix{Int64}:
 1  1  1  1
 2  2  2  2
 3  3  3  3

julia> Y = repeat(y', outer=(length(x), 1))
3×4 Matrix{Int64}:
 1  2  3  4
 1  2  3  4
 1  2  3  4

julia> X = X .+ 0.1Y
^[[B3×4 Matrix{Float64}:
 1.1  1.2  1.3  1.4
 2.1  2.2  2.3  2.4
 3.1  3.2  3.3  3.4

julia> Y = Y + 0.1cos.(X)
3×4 Matrix{Float64}:
 1.04536   2.03624  3.02675  4.017
 0.949515  1.94115  2.93337  3.92626
 0.900086  1.90017  2.90125  3.90332

julia> Z = X + cos.(Y)
3×4 Matrix{Float64}:
 1.60159  0.751184  0.306587  0.759313
 2.68208  1.83805   1.3216    1.69238
 3.72154  2.87655   2.32874   2.67636

julia> fig, ax, srf = surface(X, Y, zeros(size(Z)), color = Z, shading = NoShading, interpolate = false)

and produces something like

but it throws the following with CairoMakie v0.12.9:

Error showing value of type Makie.FigureAxisPlot:
ERROR: MethodError: no method matching iterate(::ColorTypes.RGBA{Float32})

Closest candidates are:
  iterate(::JuliaInterpreter.ExprSplitter)
   @ JuliaInterpreter ~/.julia/packages/JuliaInterpreter/7b5LK/src/construct.jl:535
  iterate(::JuliaInterpreter.ExprSplitter, ::Any)
   @ JuliaInterpreter ~/.julia/packages/JuliaInterpreter/7b5LK/src/construct.jl:535
  iterate(::Automa.Precondition)
   @ Automa ~/.julia/packages/Automa/rbg94/src/precond.jl:87
  ...

Stacktrace:
  [1] indexed_iterate(I::ColorTypes.RGBA{Float32}, i::Int64)
    @ Base ./tuple.jl:95
  [2] draw_pattern(ctx::Cairo.CairoContext, zorder::Vector{…}, shading::Bool, meshfaces::Vector{…}, ts::Vector{…}, per_face_col::CairoMakie.FaceIterator{…}, ns::Vector{…}, vs::Vector{…}, lightdir::Vec{…}, light_color::Vec{…}, shininess::Float32, diffuse::Vec{…}, ambient::Vec{…}, specular::Vec{…})
    @ CairoMakie ~/.julia/packages/CairoMakie/ramgi/src/primitives.jl:1149
  [3] draw_mesh3D(scene::Scene, screen::CairoMakie.Screen{…}, space::Symbol, transform_func::Tuple{…}, meshpoints::Vector{…}, meshfaces::Vector{…}, meshnormals::Vector{…}, per_face_col::CairoMakie.FaceIterator{…}, pos::Vec{…}, scale::Float32, rotation::StaticArraysCore.SMatrix{…}, model::StaticArraysCore.SMatrix{…}, shading::Bool, diffuse::Vec{…}, specular::Vec{…}, shininess::Float32, faceculling::Int64, clip_planes::Vector{…})
    @ CairoMakie ~/.julia/packages/CairoMakie/ramgi/src/primitives.jl:1105
  [4] draw_mesh3D(scene::Scene, screen::CairoMakie.Screen{…}, attributes::Surface{…}, mesh::GeometryBasics.Mesh{…}; pos::Vec{…}, scale::Float32, rotation::StaticArraysCore.SMatrix{…}, uv_transform::StaticArraysCore.SMatrix{…})
    @ CairoMakie ~/.julia/packages/CairoMakie/ramgi/src/primitives.jl:1015
  [5] draw_atomic(scene::Scene, screen::CairoMakie.Screen{CairoMakie.IMAGE}, primitive::Surface)
    @ CairoMakie ~/.julia/packages/CairoMakie/ramgi/src/primitives.jl:1195
  [6] draw_plot(scene::Scene, screen::CairoMakie.Screen{CairoMakie.IMAGE}, primitive::Surface{Tuple{…}})
    @ CairoMakie ~/.julia/packages/CairoMakie/ramgi/src/infrastructure.jl:129
  [7] cairo_draw(screen::CairoMakie.Screen{CairoMakie.IMAGE}, scene::Scene)
    @ CairoMakie ~/.julia/packages/CairoMakie/ramgi/src/infrastructure.jl:51
  [8] display(screen::CairoMakie.Screen{CairoMakie.IMAGE}, scene::Scene; connect::Bool)
    @ CairoMakie ~/.julia/packages/CairoMakie/ramgi/src/display.jl:43
  [9] display(screen::CairoMakie.Screen{CairoMakie.IMAGE}, scene::Scene)
    @ CairoMakie ~/.julia/packages/CairoMakie/ramgi/src/display.jl:40
 [10] display(figlike::Makie.FigureAxisPlot; backend::Module, inline::MakieCore.Automatic, update::Bool, screen_config::@Kwargs{})
    @ Makie ~/.julia/packages/Makie/8h0bl/src/display.jl:166
 [11] display(figlike::Makie.FigureAxisPlot)
    @ Makie ~/.julia/packages/Makie/8h0bl/src/display.jl:130
 [12] #invokelatest#2
    @ ./essentials.jl:892 [inlined]
 [13] invokelatest
    @ ./essentials.jl:889 [inlined]
 [14] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{…})
    @ REPL ~/.julia/juliaup/julia-1.10.4+0.x64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/REPL.jl:315
 [15] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL ~/.julia/juliaup/julia-1.10.4+0.x64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/REPL.jl:284
 [16] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL ~/.julia/juliaup/julia-1.10.4+0.x64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/REPL.jl:569
 [17] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL ~/.julia/juliaup/julia-1.10.4+0.x64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/REPL.jl:282
 [18] (::REPL.var"#do_respond#80"{…})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL ~/.julia/juliaup/julia-1.10.4+0.x64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/REPL.jl:911
 [19] #invokelatest#2
    @ ./essentials.jl:892 [inlined]
 [20] invokelatest
    @ ./essentials.jl:889 [inlined]
 [21] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit ~/.julia/juliaup/julia-1.10.4+0.x64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/LineEdit.jl:2656
 [22] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL ~/.julia/juliaup/julia-1.10.4+0.x64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/REPL.jl:1312
 [23] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ~/.julia/juliaup/julia-1.10.4+0.x64.apple.darwin14/share/julia/stdlib/v1.10/REPL/src/REPL.jl:386
Some type information was truncated. Use `show(err)` to see complete types.

(Note this is not master but I think these changes have been tagged and are part of these released versions of CairoMakie and GLMakie, right?)

[liuyxpp]

Finally, I can make the following plot with GLMakie v0.10.9:

Hoping to see this feature working in CairoMakie so that I can export PDF/SVG files.

[briochemc]

Can you share the code to produce it?

[liuyxpp]

Can you share the code to produce it?

A line of code similar to yours:

surface(X, Y, zeros(size(Z)), color = Z, shading = NoShading, interpolate = true)

X, Y, Z are generated by a simulation.

[asinghvi17]

That will already work in CairoMakie ;) , just without value interpolation (that's a vector graphics issue, nothing we can do :(

[asinghvi17]

Somehow, I couldn't get your example to work @briochemc - it fails with the same error. But this (which shouldn't be much different) works:

julia> using CairoMakie

julia> xs = 1:10
1:10

julia> ys = 1:10
1:10

julia> rp = (Makie.rotmatrix2d(π/5),) .* Point2f.(xs, ys')
10×10 Matrix{Point{2, Float64}}:
 [0.221232, 1.3968]  [-0.366554, 2.20582]  [-0.954339, 3.01484]  [-1.54212, 3.82385]   [-2.12991, 4.63287]   [-2.71769, 5.44189]   [-3.30548, 6.2509]     [-3.89327, 7.05992]   [-4.48105, 7.86894]   [-5.06884, 8.67796]
 [1.03025, 1.98459]  [0.442463, 2.7936]    [-0.145322, 3.60262]  [-0.733107, 4.41164]  [-1.32089, 5.22066]   [-1.90868, 6.02967]   [-2.49646, 6.83869]    [-3.08425, 7.64771]   [-3.67203, 8.45672]   [-4.25982, 9.26574]
 [1.83927, 2.57237]  [1.25148, 3.38139]    [0.663695, 4.19041]   [0.07591, 4.99942]    [-0.511875, 5.80844]  [-1.09966, 6.61746]   [-1.68745, 7.42647]    [-2.27523, 8.23549]   [-2.86302, 9.04451]   [-3.4508, 9.85353]
 [2.64828, 3.16016]  [2.0605, 3.96917]     [1.47271, 4.77819]    [0.884927, 5.58721]   [0.297142, 6.39623]   [-0.290644, 7.20524]  [-0.878429, 8.01426]   [-1.46621, 8.82328]   [-2.054, 9.63229]     [-2.64178, 10.4413]
 [3.4573, 3.74794]   [2.86951, 4.55696]    [2.28173, 5.36598]    [1.69394, 6.17499]    [1.10616, 6.98401]    [0.518373, 7.79303]   [-0.0694118, 8.60205]  [-0.657197, 9.41106]  [-1.24498, 10.2201]   [-1.83277, 11.0291]
 [4.26632, 4.33573]  [3.67853, 5.14475]    [3.09075, 5.95376]    [2.50296, 6.76278]    [1.91518, 7.5718]     [1.32739, 8.38081]    [0.739605, 9.18983]    [0.15182, 9.99885]    [-0.435965, 10.8079]  [-1.02375, 11.6169]
 [5.07533, 4.92351]  [4.48755, 5.73253]    [3.89976, 6.54155]    [3.31198, 7.35056]    [2.72419, 8.15958]    [2.13641, 8.9686]     [1.54862, 9.77762]     [0.960837, 10.5866]   [0.373052, 11.3956]   [-0.214734, 12.2047]
 [5.88435, 5.5113]   [5.29657, 6.32032]    [4.70878, 7.12933]    [4.12099, 7.93835]    [3.53321, 8.74737]    [2.94542, 9.55638]    [2.35764, 10.3654]     [1.76985, 11.1744]    [1.18207, 11.9834]    [0.594283, 12.7925]
 [6.69337, 6.09908]  [6.10558, 6.9081]     [5.5178, 7.71712]     [4.93001, 8.52614]    [4.34223, 9.33515]    [3.75444, 10.1442]    [3.16666, 10.9532]     [2.57887, 11.7622]    [1.99109, 12.5712]    [1.4033, 13.3802]
 [7.50238, 6.68687]  [6.9146, 7.49589]     [6.32681, 8.3049]     [5.73903, 9.11392]    [5.15124, 9.92294]    [4.56346, 10.732]     [3.97567, 11.541]      [3.38789, 12.35]      [2.8001, 13.159]      [2.21232, 13.968]

julia> surface(first.(rp), last.(rp), zeros(length(xs), length(ys)); color = Makie.peaks(length(xs)), shading = NoShading)