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:
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.
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: :wink:](https://emoji.discourse-cdn.com/twitter/wink.png?v=12)
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
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)
fig
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: :slight_smile:](https://emoji.discourse-cdn.com/twitter/slight_smile.png?v=12)