Построение звездных карт в экваториальной системе координат

Я не могу повторить ваш код, чтобы показать, где именно вам нужны исправления. Matplotlib позволяет использовать полярные координаты, для этого вам нужно применить оси и сетки, как в этих документах: matplotlib docs

Если вы можете упростить код, не используя файлы (только список чисел или строк, я не знаю), я могу попытаться это исправить.

Ниже приведен простой пример того, как это можно сделать:

    #!/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
import numpy as np

from matplotlib.projections import PolarAxes
from matplotlib.transforms import Affine2D
from mpl_toolkits.axisartist import Axes, HostAxes, angle_helper
from mpl_toolkits.axisartist.grid_helper_curvelinear import \
    GridHelperCurveLinear

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

# Query object from Simbad
OBJECT = "Polaris"
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"]

# <- Start point code implementetion from docs
tr = Affine2D().scale(np.pi/180, 1) + PolarAxes.PolarTransform()
grid_locator1 = angle_helper.LocatorDMS(12)

grid_helper = GridHelperCurveLinear(
    tr, grid_locator1=grid_locator1)

ax = fig.add_subplot(
    1, 1, 1, axes_class=HostAxes, grid_helper=grid_helper)
ax.axis["right"].major_ticklabels.set_visible(True)
ax.axis["top"].major_ticklabels.set_visible(True)
ax.axis["right"].get_helper().nth_coord_ticks = 0
ax.axis["bottom"].get_helper().nth_coord_ticks = 1

ax.set_aspect(1)
ax.grid(True, zorder=0)
ax.grid(True, color='white', linestyle='dotted')
# <- End point code implementetion from docs

# 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
# Also some changes below
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_aspect('equal')

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

# Set the coordinate grid
ax.set_xlabel('RA (hours)')
ax.set_ylabel('Dec (degrees)')
ax.xaxis.set_major_formatter('hh:mm:ss')
ax.yaxis.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')

Может я что-то пропустил, но вроде ок

Я мог найти решение на другом форуме. Итак, я опубликую ответ. Это было проще, чем кажется!

Существенная проблема заключалась в моих координатных шкалах WCS, они были определены неправильно. На самом деле, при использовании поля зрения 30 изображение не имело поля обзора 30 градусов, оно было меньше. Этот момент помог найти ответ.

Возникла проблема с используемыми единицами измерения. Хотя я рассчитывал все в радианах, WCSAxes работает в пикселях, поэтому мне нужно было определить CDELT равным 1 радиану на пиксель.

Итак, исправленная строка:

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

Это окончательное изображение:

Я очень-очень этому рад!