import astropy.units as u
from astropy.coordinates import SkyCoord
from gala.coordinates import GreatCircleICRSFrame
import matplotlib.pyplot as plt

pt1 = SkyCoord(ra=170.*u.deg, dec=23.18*u.deg)
pt2 = SkyCoord(ra=125.7*u.deg, dec=-72.2*u.deg)
gcf2 = GreatCircleICRSFrame.from_endpoints(pt1, pt2)

grid_c = SkyCoord(phi1=np.linspace(0, 360, 128)*u.deg, phi2=0*u.deg,
                 frame=gcf2)
grid_c = grid_c.galactic

plt.figure(figsize=(6, 4))
plt.plot(grid_c.l.degree[grid_c.l.degree.argsort()],
        grid_c.b.degree[grid_c.l.degree.argsort()],
        ls='-', marker='')
plt.xlabel('$l$ [deg]')
plt.ylabel('$b$ [deg]')
plt.tight_layout()