如何解决对于给定的纬度,如何在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()
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。