I am new to Skyfield and have been trying to automate the creation of full-sky star charts, given a location on earth and a time. I have attempted to create a point that corresponds to a given location's zenith and then observe that point, as in the comet example found in the documentation: https://rhodesmill.org/skyfield/example-plots.html#drawing-a-finder-chart-for-comet-neowise
However, the method I've used generates an ICRF position, which is apparently incompatible with the observe() method. I've been reading about positions and reference frames and haven't quite wrapped my head around how they're handled in Skyfield. Is there an easy way to convert an ICRF position to make it compatible with observe(), or is there an altogether simpler way to generate a star chart that centers on the zenith?
Here's my code snippet, largely taken from the example at the link:
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.collections import LineCollection
from skyfield.api import Star, load, N, W, wgs84, Topos
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from skyfield.data import hipparcos, mpc, stellarium
from skyfield.projections import build_stereographic_projection
from skyfield.positionlib import Astrometric, position_of_radec
from skyfield.vectorlib import VectorFunction, VectorSum
ts = load.timescale()
t = ts.utc(2019,9,13,20)
lat = 41.0
lon = 111.0
geographic = wgs84.latlon(latitude_degrees=lat*N, longitude_degrees=lon*W)
observer = geographic.at(t)
zenith = observer.from_altaz(alt_degrees=90, az_degrees=0)
zenith_ra, zenith_dec, dist = zenith.radec()
eph = load('de421.bsp')
sun = eph['sun']
earth = eph['earth']
viewloc = earth + wgs84.latlon(lat*N, lon*W)
# The Hipparcos mission provides our star catalog.
with load.open(hipparcos.URL) as f:
stars = hipparcos.load_dataframe(f)
url = ('https://raw.githubusercontent.com/Stellarium/stellarium/master'
'/skycultures/western_SnT/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]
# We will center the chart on the zenith.
barycentric = earth.at(t)
center_position = position_of_radec(zenith_ra.hours, zenith_dec.degrees)
astrometric = Astrometric(center_position.position.au) # Not sure what is needed as an argument here
center = barycentric.observe(astrometric)
projection = build_stereographic_projection(center)
field_of_view_degrees = 45.0
limiting_magnitude = 7.0
Try using the Zenith position directly as the center of your plot:
projection = build_stereographic_projection(zenith)
See whether that puts the right stars at the center of the resulting diagram. If not, try pasting in the result as an image, and we can explore further!