对于给定的纬度，如何在python中用赤道坐标绘制天球的正投影？

问题描述

我试图获得天球的正投影，从某个纬度看到的赤道坐标，如下图所示：

（从 Skychart/Cartes du ciel 获得的网格）

这张图片是 Skychart/Cartes du ciel 的打印件，显示了观察者在 23°S 纬度的赤道网格。我希望能够在 Python 中重现完全相同的图像（除了深蓝色背景）。我的第一次尝试是使用cartopy，设置中心纬度为-23，如下：

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
ax = plt.axes(projection=ccrs.Orthographic(central_latitude=-23))
ax.gridlines()
plt.show()

但生成的图片是这样的：

从南极的位置来看，我认为在 cartopy 中将中心纬度设置为观察者的纬度并不能解决我的问题。是否有另一种方法，无论是使用 cartopy 还是其他包（也许是 AstroPy？ - 我从未使用过它）来获得与 Python 中的 Skychart（图 1）相同的图？

解决方法

首先，您的第一个图像是方位角等距投影。因此，它与您的第二个图（正交投影）完全不同。要使用 Cartopy 获得像这样的图（第一张图像），需要遵循一些有趣的步骤。下面是带有注释的代码，它生成了我认为效果不错的输出图。

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.path import Path

import matplotlib.path as mpath
import numpy as np

r_limit = 20037508 #from: ax.get_ylim() of full extent

# this makes circle for clipping the plot
pts = [] #unit circle vertices
cds = [] #path codes
numps = 32
for ix,ea in enumerate(np.linspace(0,2*np.pi,numps)):
    #print(ea)
    xi = np.cos(ea)
    yi = np.sin(ea)
    pts.append([xi,yi])
    if (ix==0):
        # start
        cds.append(1)
    elif (ix==numps-1):
        # close
        cds.append(79)
    else:
        cds.append(4)

# make them np.array for easy uses
vertices = np.array(pts) 
codes = np.array(cds)

# manipulate them to create a required clip_path
scale = r_limit*0.5
big_ccl = mpath.Path(vertices*scale,codes)
clippat = plt.Polygon(big_ccl.vertices[:,:],visible=True,fill=False,ec='red')

# create axis to plot `AzimuthalEquidistant` projection
# this uses specific `central_latitude`
ax = plt.axes(projection=ccrs.AzimuthalEquidistant(central_latitude=-23))

# add the clip_path
ax.add_patch(clippat)

# draw graticule (of meridian and parallel lines)
# applying clip_path to get only required extents plotted
ax.gridlines(draw_labels=False,crs=ccrs.PlateCarree(),xlocs=range(-180,180,30),ylocs=range(-80,90,20),clip_path=clippat)

# specify radial extents,use them to set limits of plot
r_extent = r_limit/(2-0.05)  # special to the question
ax.set_xlim(-r_extent,r_extent)
ax.set_ylim(-r_extent,r_extent)

ax.set_frame_on(False)  #hide the rectangle frame

plt.show()