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

Question

I am trying to obtain an ortographic projection of the celestial sphere, with equatorial coordinates, as seen from a certain latitude, as in the following picture:

(Grid obtained from Skychart/Cartes du ciel)

This image is a print of Skychart/Cartes du ciel, showing the equatorial grid for an observer at 23°S latitude. I want to be able to reproduce the exact same image in Python (apart from the dark blue background). My first attempt was to use CartoPy, setting the central latitude as -23, as follows:

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

but the resulting picture looks like this:

From the position of the south pole, I believe setting the central latitude to the observer's latitude in CartoPy does not solve my problem. Is there another way, either with CartoPy or another package (maybe AstroPy? - I have never used it) to obtain the same plot as Skychart (Image 1) in python?

question from:https://stackoverflow.com/questions/66048195/how-to-plot-an-ortographic-projection-of-the-celestial-sphere-with-equatorial-co

Answer 1

First of all, your first image is Azimuthal Equidistant Projection. So that, it is quite different from your second plot (Orthographic projection). To get the plot (first image) like that using Cartopy requires some steps that are interesting to follow. Here is the code with comments that produces the output plot that I consider a good result.

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