Plotting HEALPix map onto mollweide projection with GeoMakie

Hi, I am trying to create a figure like this:

I’m using this all-sky map: LAMBDA - Foreground - 2014 Reprocessed Haslam Map Downloads

This is my best attempt at the moment

using GeoMakie, GLMakie
using Healpix

function main()
    f = Figure(size=(800, 500))
    ax = GeoAxis(f[1, 1], xticks=-180:30:180; dest="+proj=moll")

    ylims!(ax, -90, 90)
    xlims!(ax, 180, -180)
    ax.xreversed = true

    healpix_map = readMapFromFITS("./haslam408_ds_Remazeilles2014.fits", 1, Float32)
    moll = mollweide(healpix_map)

    meshimage!(ax, -180 .. 180, -90 .. 90, moll[1]; npoints=1000)
    arc!(ax, (0, -60), 10, -pi, pi; linewidth=5)
    arc!(ax, (-150, -60), 10, -pi, pi; linewidth=5)


Which gets me this:

I have tried including source="+proj=moll" inside the meshimage!, but the image just disappears from the axis completely. How should I go about fixing the white space? I would like to keep it on a GeoAxis so I can easily include more arcs! at different positions and diff radius, like shown in the code.


I have seen issues with the meshimage recipe and transformations before, am looking into it. In the meanwhile, does surface work instead of meshimage?

Also, could you give me an MWE, maybe with generated random data that I can run on my machine?

Also, if you set source="+proj=moll", then your bounds also have to be in the space of the Mollweide projection :wink:

So the problem you are facing is not that there is something wrong with the projection, but rather that you are trying to project an image that has already been projected. For example, with a minimal example from the Healpix.jl documentation,


will yield

Now, plotting with meshimage will cause the entire contents of that axis, along with all the blank space, to be translated to the Mollweide projection again.

What I suspect you really want is:

rectimg = equirectangular(healpix_map)
meshimage!(ax, -180 .. 180, -90 .. 90, reverse(rectimg[1]; dims = 1); npoints=1000)

which will yield (with an extra lines!(ax, GeoMakie.coastlines(); color = :black)):

1 Like

See also Add an example of how to use Healpix.jl with GeoMakie. by asinghvi17 · Pull Request #246 · MakieOrg/GeoMakie.jl · GitHub, which adds an example of how to use Healpix with GeoMakie.

In that context, an adapted version of your code is:

healpix_map = readMapFromFITS("./haslam408_ds_Remazeilles2014.fits", 1, Float32)
erect, _, _ = equirectangular(healpix_map)

fig = Figure(size=(800, 500))
ax = GeoAxis(fig[1, 1], xticks=-180:30:180; dest="+proj=moll")

ylims!(ax, -90, 90)
xlims!(ax, -180, 180)

meshimage!(ax, -180 .. 180, -90 .. 90, reverse(erect; dims = 1); colorrange = Makie.PlotUtils.zscale(erect), npoints=720)
arc!(ax, (0, -60), 10, -pi, pi; linewidth=5)
arc!(ax, (-150, -60), 10, -pi, pi; linewidth=5)

Changes are:

  • Use equirectangular instead of mollweide as a projection
  • Use npoints=720, since that’s the number of pixels that Healpix’s projection functions return. Using more isn’t really useful, and one can get away with even less (~300) points in meshimage if using GLMakie.
  • Colorrange changed to use zscale (not sure if you want to use this or not)

Thank you, this is exactly what I wanted! I’ll just do some scaling of the data now and add some padding on the ticks. Honestly tried to do this in python as well, but didn’t even manage to get an image. Hopefully I can use Julia and Makie for plots more in the future!

FYI @AtomsForHire, I just edited my last post in this thread with the results of that code using the corrected meshimage recipe - the scaling looks a lot better now!

1 Like

Yep looks much nicer now! Thanks for being so active in here as well :slight_smile: