Plotting HEALPix map onto mollweide projection with GeoMakie

Specific Domains Visualization
, ,
1

Hi, I am trying to create a figure like this: https://www.researchgate.net/publication/353085441/figure/fig1/AS:1043304233848839@1625754521890/All-sky-map-showing-the-best-fit-positions-and-90-containment-regions-approximated-to.png

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)
    wait(display(f))
end

main()

Which gets me this:

image
image1920×1197 70.3 KB

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.

Thanks.

2

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:

3

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,

heatmap(moll[1])

will yield

download
download1200×900 45.9 KB

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)):

download-1
download-11200×900 259 KB

1 Like
4

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)
fig

iTerm2.O7IJrv
iTerm2.O7IJrv1600×1000 87.4 KB

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)
5

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!

6

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
7

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