Plotting star maps with equatorial coordinates system

I'm trying to generate star maps with the equatorial coordinates system (RAJ2000 and DEJ2000). However, I only get a grid system where meridians and parallels are in parallel, while parallels should be curved and meridians should converge to the north celestial pole and the south ceestial pole.

I'm using some Python modules: matplotlib, skyfield (for the stereographic projection), astroquery (so I can target any object in the deep space) and astropy.

This is my code:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Generate a skymap with equatorial grid"""

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.collections import LineCollection
from skyfield.api import Star, load
from skyfield.data import hipparcos, stellarium
from skyfield.projections import build_stereographic_projection
from astroquery.simbad import Simbad
from astropy.coordinates import SkyCoord
import astropy.units as u
from astropy.wcs import WCS
from astropy.visualization.wcsaxes import WCSAxes

# Design
plt.style.use("dark_background")
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman']

# Query object from Simbad
OBJECT = "Alioth"
FOV = 30.0
MAG = 6.5

TABLE = Simbad.query_object(OBJECT)
RA = TABLE['RA'][0]
DEC = TABLE['DEC'][0]
COORD = SkyCoord(f"{RA} {DEC}", unit=(u.hourangle, u.deg), frame='fk5')

print("RA is", RA)
print("DEC is", DEC)

ts = load.timescale()
t = ts.now()

# An ephemeris from the JPL provides Sun and Earth positions.
eph = load('de421.bsp')
earth = eph['earth']

# Load constellation outlines from Stellarium
url = ('https://raw.githubusercontent.com/Stellarium/stellarium/master'
       '/skycultures/modern_st/constellationship.fab')

with load.open(url) as f:
    constellations = stellarium.parse_constellations(f)

edges = [edge for name, edges in constellations for edge in edges]
edges_star1 = [star1 for star1, star2 in edges]
edges_star2 = [star2 for star1, star2 in edges]

# The Hipparcos mission provides our star catalog.
with load.open(hipparcos.URL) as f:
    stars = hipparcos.load_dataframe(f)

# Center the chart on the specified object's position.
center = earth.at(t).observe(Star(ra_hours=COORD.ra.hour, dec_degrees=COORD.dec.degree))
projection = build_stereographic_projection(center)

# Compute the x and y coordinates that each star will have on the plot.
star_positions = earth.at(t).observe(Star.from_dataframe(stars))
stars['x'], stars['y'] = projection(star_positions)

# Create a True/False mask marking the stars bright enough to be included in our plot.
bright_stars = (stars.magnitude <= MAG)
magnitude = stars['magnitude'][bright_stars]
marker_size = (0.5 + MAG - magnitude) ** 2.0

# The constellation lines will each begin at the x,y of one star and end at the x,y of another.
xy1 = stars[['x', 'y']].loc[edges_star1].values
xy2 = stars[['x', 'y']].loc[edges_star2].values
lines_xy = np.rollaxis(np.array([xy1, xy2]), 1)

# Define the limit for the plotting area
angle = np.deg2rad(FOV / 2.0)
limit = np.tan(angle)  # Calculate limit based on the field of view

# Build the figure with WCS axes
fig = plt.figure(figsize=[6, 6])
wcs = WCS(naxis=2)
wcs.wcs.crpix = [1, 1]
wcs.wcs.cdelt = np.array([-FOV / 360, FOV / 360])
wcs.wcs.crval = [COORD.ra.deg, COORD.dec.deg]
wcs.wcs.ctype = ["RA---STG", "DEC--STG"]

ax = fig.add_subplot(111, projection=wcs)

# Draw the constellation lines
ax.add_collection(LineCollection(lines_xy, colors='#ff7f2a', linewidths=1, linestyle='-'))

# Draw the stars
ax.scatter(stars['x'][bright_stars], stars['y'][bright_stars],
           s=marker_size, color='white', zorder=2)

ax.scatter(RA, DEC, marker='*', color='red', zorder=3)

angle = np.pi - FOV / 360.0 * np.pi
limit = np.sin(angle) / (1.0 - np.cos(angle))

# Set plot limits
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_aspect('equal')

# Add RA/Dec grid lines
ax.coords.grid(True, color='white', linestyle='dotted')

# Set the coordinate grid
ax.coords[0].set_axislabel('RA (hours)')
ax.coords[1].set_axislabel('Dec (degrees)')
ax.coords[0].set_major_formatter('hh:mm:ss')
ax.coords[1].set_major_formatter('dd:mm:ss')

# Title
ax.set_title(f'Sky map centered on {OBJECT}', color='white', y=1.04)

# Save the image
FILE = "chart.png"
plt.savefig(FILE, dpi=100, facecolor='#1a1a1a')

And this is the resulting image:

As you can see, the grid (parallels and meridians) are totally parallel. However, my goal is to achieve this grid:

In this case, I got the right WCS from a FITS image in the DSS survey. It's automatic. However, for plotting star maps I need to create a simulation of that, working fine with the labels and the coordinates system, not as a background image or something else.

I could find the solution in another forum. So I'm gonna post the answer. It was easier than it seems!

The essential problem was in my WCS coordinate scales, they were defined incorrectly. Actually, using the FOV 30, the image didn't have a FOV of 30 degrees, it was smaller. This point helped to find the answer.

There was a problem in the units used. While I was plltting everything in radians, WCSAxes works in pixels, so I needed to define CDELT to be 1 radian per pixel.

So the line corrected is this one:

wcs.wcs.cdelt = np.array([-360 / np.pi, 360 / np.pi])

This is the final image:

I'm very very happy with this!