I'm attempting to draw a scatter plot of astronomic sources projected on a full sky map, using two coordinate grids: galactic and equatorial. Generating a single grid plot (using one of the pre-defined projections) is easy:
import matplotlib.pyplot as plt
import numpy as np
# Random galactic coordinates
l_rad = np.random.uniform(-np.pi, np.pi, 100)
b_rad = np.random.uniform(-np.pi / 2., np.pi / 2., 100)
plt.figure()
plt.subplot(111, projection="aitoff")
plt.grid(True)
plt.plot(ra_rad, dec_rad, 'o', markersize=2, alpha=0.3)
plt.show()
What I need though is a map with two grids superimposed: one for galactic coordinates (which is already there) and another one for equatorial. That is, a second grid like the one in the image below (notice the hour convention instead of decimal degrees for the right ascension) should be superimposed onto the first grid:
I think the Kapteyn package can do something similar (image below) but I'd prefer to avoid it if possible (I find it rather complicated to use).