I have three healpy maps (corresponding to three Stokes parameters at different frequencies) and I am trying to plot them together in the same map, to highlight differences. I want to do something like what is done in the big map in the first plot: https://sci.esa.int/web/planck/-/56323-planck-legacy-archive-foreground-component-maps . Is there a way of doing this? I have appreciated that Healpy does not have any functionality to create RGB maps from it.
You cannot do this in healpy, but you can probably do with matplotlib.
You can project a map using a mollweide projection:
https://matplotlib.org/stable/gallery/subplots_axes_and_figures/geo_demo.html
So you can get the 3 maps, normalize them appropriately (not sure if between 0 and 1 or 0 and 256), then combine them in a single 2D array where first row is R, second G and last B.
Then you can plot it with imshow, see that the input of imshow can be a RGB array:
https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.imshow.html