R: plotting sky map in Mollweide projection in galactic coordinate system

Question

I need to produce the star map with constellations etc in the Mollweide projection (elliptical projection giving 360 view angle, used in plotting night sky). I found a recipe at https://kimnewzealand.github.io/2019/02/21/celestial-maps/ with the use of sf package and converting the default EPSG:4326 data of the sky objects into Mollweide projection.

At some stage the data is converted to the Mollweide projection using the command:

constellation_lines_sf_trans<- st_transform(constellation_lines_sf_trans, crs = "+proj=moll")

The resulting image, reproduced along the lines as described in the link, looks like:

It is fine, however, the coordinate system is equatorial, that is basically with the same rotation axis as all coordinate systems on the Earth, like WGS84 (North Pole upwards). For example, the Milky Way is shown on this plot, going at some angle 60 degrees. We need the so called galactic coordinates: this is the coordinate plane coinciding with the plane of our Galaxy. So, Milky way here would be just a horizontal line of the ellipse axis. For example, the solution found elsewhere, seems to use the same technique, but the code is not given there:

Here Milky Way is a horizontal line, and the North Pole is in the upper left corner (denoted as np; for example, here one can see the distorted recognizable constellations of Ursa Major/Minor around the North Pole). I would take this image, but there is a blind spot (showing the blind zone of an observatory which cannot reach this region in the sky), so I would like to reproduce this image: constellations + Mollweide projection + "galactic" orientation of the reference frame.

We are able to convert between variety of coordinate systems in R packages. It seems that most GIS tools use various flavors of Earth-related coordinate systems and projections, based to the rotation of the Earth (North Pole up), for majority of applications, needed for GIS. The question is whether it is possible to load and convert to a predefined galactic coordinate system (or, for example, to the ecliptic system), or to perform this conversion on the fly in the scripts with manual conversion of star data

Answer 1

EDIT: Actually, after further research, it appears that it all comes down to having your projection do a rotation. This also happens in the code on the interactive example I mention farther below. To make it center of the Galactic center, you need to rotate by [93.5949, 28.9362, -58.5988], which specify the [lambda, phi, gamma] rotation angles in degrees about each spherical axis. (when you rotate, there's no need to convert your ra or dec coordinates to galactic coordinates anymore.)

I don't know enough about mapping in R to be able to say if it's possible to specify a rotation on a projection, but here's an amazing example that shows this process with d3.js (and for those truly interested, shows where the angles come from). In case it's not possible, perhaps the route below would still be viable.

---

I was investigating the same thing, and I might've found something. I think that you first need to convert your equatorial ra (right ascension) and dec (declination) coordinates to galactic coordinates. And then apply a (Mollweide) projection. I'm not sure if that is completely correct since my case was slightly different, but at least this worked for me:

• The dataset has rows of ra &dec in equatorial coordinates (given in degrees in my case)
• Using the euler function of the astrolibR package, I calculate the Galactic longitude gl and latitude gb (ps: you can also use the glactc function):

data$gl <- euler(data$ra, data$dec, select=1)$ao
data$gb <- euler(data$ra, data$dec, select=1)$bo

• Next, using these new coordinates I apply the Aitoff projection that is also part of the astrolibR package to get back x and y coordinates:

data$x <- -aitoff(data$gl, data$gb)$x
data$y <- aitoff(data$gl, data$gb)$y

which I can then plot

ggplot(data, aes(x, y)) + geom_point(shape=16, size = 0.1, alpha = 0.2) + coord_fixed()

(the image below is based on my own dataset of observations, the "darker line" follows the ecliptic line, and the two blobs in the lower right are the Large & Small Magellanic Clouds) I found it useful to compare against this interactive map and setting the coordinates to "galactic", centering on 0,0 and then trying different kinds of projections.

Perhaps you can try applying the Mollweide projection with the gl and gb coordinates instead?