# Plotting HEALPix map onto mollweide projection with GeoMakie

**URL:** https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805
**Category:** Visualization
**Tags:** question, plotting, makie
**Created:** [June 18, 2024, 12:48pm UTC](https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805 "2024-06-18T12:48:30Z")
**Posts on this page:** 7
**Page:** 1

<div class="post-metadata">

### Author: ![AtomsForHire](https://avatars.discourse-cdn.com/v4/letter/a/278dde/32.png) [@AtomsForHire](https://discourse.julialang.org/u/AtomsForHire)
#### Post date: [June 18, 2024, 12:48pm UTC](https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805/1 "2024-06-18T12:48:30Z")

</div>

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](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](https://lambda.gsfc.nasa.gov/product/foreground/fg_2014_haslam_408_get.html)

This is my best attempt at the moment

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/5/9/59b702b87b2cdfcc36ff1a09b9a5a1c73950fe34.jpeg)

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.

---

<div class="post-metadata">

### Author: ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)
#### Post date: [June 18, 2024, 6:41pm UTC](https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805/2 "2024-06-18T18:41:52Z")

</div>

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 😉

---

<div class="post-metadata">

### Author: ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)
#### Post date: [June 18, 2024, 8:28pm UTC](https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805/3 "2024-06-18T20:28:51Z")

</div>

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,

```julia
heatmap(moll[1])

```

will yield

 ![download](https://global.discourse-cdn.com/julialang/original/3X/5/3/531e17cdd5216b905075cbfc5180a082701b9d9a.png)  
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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/2/2/2206da1c5fc103b78c6a1077621af367337763df.png)

---

<div class="post-metadata">

### Author: ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)
#### Post date: [June 18, 2024, 8:31pm UTC](https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805/4 "2024-06-18T20:31:03Z")

</div>

See also [Add an example of how to use Healpix.jl with GeoMakie. by asinghvi17 · Pull Request #246 · MakieOrg/GeoMakie.jl · GitHub](https://github.com/MakieOrg/GeoMakie.jl/pull/246), which adds an example of how to use Healpix with GeoMakie.

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/d/7/d7c7060fe3ed6cfd45d93c9bc8ea222b34a29055.jpeg)

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)

---

<div class="post-metadata">

### Author: ![AtomsForHire](https://avatars.discourse-cdn.com/v4/letter/a/278dde/32.png) [@AtomsForHire](https://discourse.julialang.org/u/AtomsForHire)
#### Post date: [June 19, 2024, 2:45am UTC](https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805/5 "2024-06-19T02:45:06Z")

</div>

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!

---

<div class="post-metadata">

### Author: ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)
#### Post date: [June 21, 2024, 12:24am UTC](https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805/6 "2024-06-21T00:24:45Z")

</div>

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!

---

<div class="post-metadata">

### Author: ![AtomsForHire](https://avatars.discourse-cdn.com/v4/letter/a/278dde/32.png) [@AtomsForHire](https://discourse.julialang.org/u/AtomsForHire)
#### Post date: [June 22, 2024, 12:59am UTC](https://discourse.julialang.org/t/plotting-healpix-map-onto-mollweide-projection-with-geomakie/115805/7 "2024-06-22T00:59:19Z")

</div>

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