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

Question

我正在尝试获取天球的正射投影，具有赤道坐标，从某个纬度看，如下图所示：

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

此图片是 Skychart/Cartes du ciel 的印刷品，显示了南纬 23° 处观察者的赤道网格。我希望能够在 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）相同的图？

Answer 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()

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

How to plot an ortographic projection of the celestial sphere with equatorial coordinates in python, for a given latitude?

python

astronomy

cartopy

astropy