PGFPlots Healpix / Mollweide / Aitoff projections?

I really like PGFPlotsX.jl to create plots which seamlessly integrate into papers but I now I need to make a bunch of Healpix Mollweide projections to visualise data spread on the surface on a sphere (essentially skymaps) and I have not really figured out how. I searched the web but it seems that either nobody is talking about, or nobody is using PGFPlots to do this :wink:

The package Healpix.jl already provides a recipe for Plots.jl which works:

julia> using Healpix

julia> using Plots

julia> m = HealpixMap{Float64, RingOrder}(32);

julia> m[ang2pixRing(m.resolution, π/4, 0.0)] = 250
250

julia> m[ang2pixRing(m.resolution, 0.0, 0.0)] = 200
200

julia> m[ang2pixRing(m.resolution, 0.42, 0.23)] = 300
300

julia> plot(m)

and produces this beauty:

But I need to adjust a lot of things and annotate stuff etc. which I find much more easier (and possible) in PGFPlots via PGFPlotsX.jl. I played around with Plots.jl but I find myself struggling really with finding the right point to override defaults and modify specific aspects of the image.

So my question: does anybody here by chance has some experience in creating healpix/mollweide/aitoff projections with PGFPlots?

I fear the answer is “no, take what’s there (Plots.jl) and try to improve”, but I’d give up so much freedom (and free time :laughing: ) by not using PGF.

Edit: let me add a nicer example output. The real fun begins with drawing lines, circles etc., modifying the colorbar and so on. Maybe it’s perfectly doable in Plots.jl, but I am a bit lost:

OK, for the sake of completeness, I manage to get something but I’ll stop working on it since I think that the output format will not be suitable for a nice user experience :wink: I thought that some renderers would somehow manage to simplify and pixelise the output but apparently everything is vectorised, which means that the output file is huge (in my case 10MB) and on slow system the scrolling is causing a lot of stress on the hardware.

Anyways, here is the rough idea:

using PGFPlotsX

out = rand(100, 100)  # this is basically a heatmap-like structure

fig = @pgf Axis(
    {
        view = (0, 90),
        colorbar,
        "colormap/viridis"
    },
    Plot3(
        {
            surf,
            shader = "flat",
        },
        Table(TableData(1:size(out)[1], 1:size(out)[2], out)))
)

Screenshot 2024-03-12 at 19.25.36

To display a HealpixMap from Healpix.jl, we need to do:

m = HealpixMap{Float64, RingOrder}(32);

# fill the map
# ...

img, mask, anymasked = mollweide(m, Dict())

and img is our out above.

And voila, the skymap from above: