Note
Go to the end to download the full example code.
Star Globe Explorer — a linked celestial sphere and sky map#
A celestial globe wired to a flat star map, with a sidereal clock.
Left — the celestial sphere. Whatever you turn toward the camera is your zenith: the point directly overhead.
Right — the same sky as an equirectangular map, with the zenith marked and the horizon drawn as the great circle 90° away from it. Everything inside the horizon is above your head right now; everything outside is below the ground.
Bottom — a sidereal clock and a latitude control.
The coupling runs both ways. Orbit the globe and the clock and latitude follow, because a zenith is a time and a place: the right ascension overhead equals the local sidereal time, and its declination equals your latitude. Drag either control and the globe turns to match.
Star positions are real (approximate J2000, rounded to about an arcminute), carried over from Textured 3-D surfaces — a celestial globe; the Milky Way band is placed by true galactic latitude.
Note
The globe is the sky seen from outside, the way a physical celestial globe is built — so the constellations read mirrored compared with lying on your back and looking up. The flat map is the conventional inside-out view.
import numpy as np
import anyplotlib as apl
# ── Star catalogue ───────────────────────────────────────────────────────────
# name: (right ascension in hours, declination in degrees, visual magnitude).
# Swapping in a full catalogue (e.g. HYG, ~9k stars to naked-eye magnitude)
# is a change to this table alone — everything below is generic.
STARS = {
"Sirius": (6.752, -16.72, -1.46), "Canopus": (6.399, -52.70, -0.74),
"Rigil Kent": (14.660, -60.83, -0.27), "Arcturus": (14.261, 19.18, -0.05),
"Vega": (18.615, 38.78, 0.03), "Capella": (5.278, 46.00, 0.08),
"Rigel": (5.242, -8.20, 0.13), "Procyon": (7.655, 5.22, 0.34),
"Achernar": (1.629, -57.24, 0.46), "Betelgeuse": (5.919, 7.41, 0.50),
"Hadar": (14.064, -60.37, 0.61), "Altair": (19.846, 8.87, 0.77),
"Acrux": (12.443, -63.10, 0.77), "Aldebaran": (4.599, 16.51, 0.85),
"Spica": (13.420, -11.16, 1.04), "Antares": (16.490, -26.43, 1.09),
"Pollux": (7.755, 28.03, 1.14), "Fomalhaut": (22.961, -29.62, 1.16),
"Deneb": (20.690, 45.28, 1.25), "Mimosa": (12.795, -59.69, 1.25),
"Regulus": (10.140, 11.97, 1.35), "Adhara": (6.977, -28.97, 1.50),
"Castor": (7.577, 31.89, 1.58), "Shaula": (17.560, -37.10, 1.62),
"Gacrux": (12.520, -57.11, 1.63), "Bellatrix": (5.418, 6.35, 1.64),
"Elnath": (5.438, 28.61, 1.65), "Miaplacidus": (9.220, -69.72, 1.67),
"Alnilam": (5.604, -1.20, 1.69), "Alnair": (22.137, -46.96, 1.74),
"Alioth": (12.900, 55.96, 1.76), "Alnitak": (5.679, -1.94, 1.77),
"Dubhe": (11.062, 61.75, 1.79), "Mirfak": (3.405, 49.86, 1.79),
"Wezen": (7.140, -26.39, 1.83), "Kaus Aust.": (18.403, -34.38, 1.85),
"Alkaid": (13.792, 49.31, 1.85), "Avior": (8.375, -59.51, 1.86),
"Sargas": (17.622, -43.00, 1.86), "Menkalinan": (5.992, 44.95, 1.90),
"Atria": (16.811, -69.03, 1.91), "Alhena": (6.628, 16.40, 1.93),
"Peacock": (20.427, -56.74, 1.94), "Polaris": (2.530, 89.26, 1.98),
"Mirzam": (6.378, -17.96, 1.98), "Alphard": (9.460, -8.66, 2.00),
"Hamal": (2.120, 23.46, 2.00), "Diphda": (0.727, -17.99, 2.04),
"Nunki": (18.921, -26.30, 2.05), "Menkent": (14.112, -36.37, 2.06),
"Mirach": (1.162, 35.62, 2.06), "Alpheratz": (0.140, 29.09, 2.06),
"Rasalhague": (17.582, 12.56, 2.08), "Kochab": (14.845, 74.16, 2.08),
"Algieba": (10.333, 19.84, 2.08), "Saiph": (5.796, -9.67, 2.09),
"Tiaki": (22.711, -46.88, 2.11), "Algol": (3.136, 40.96, 2.12),
"Denebola": (11.818, 14.57, 2.14), "Muhlifain": (12.692, -48.96, 2.20),
"Aspidiske": (9.285, -59.28, 2.21), "Alphecca": (15.578, 26.71, 2.22),
"Sadr": (20.370, 40.26, 2.23), "Mizar": (13.399, 54.93, 2.23),
"Eltanin": (17.943, 51.49, 2.23), "Suhail": (9.133, -43.43, 2.23),
"Schedar": (0.675, 56.54, 2.24), "Mintaka": (5.533, -0.30, 2.25),
"Caph": (0.153, 59.15, 2.28), "Gamma Cas": (0.945, 60.72, 2.47),
"Merak": (11.030, 56.38, 2.37), "Enif": (21.737, 9.88, 2.38),
"Phecda": (11.897, 53.69, 2.44), "Ruchbah": (1.430, 60.24, 2.68),
"Delta Cru": (12.252, -58.75, 2.79), "Alcyone": (3.792, 24.11, 2.87),
"Megrez": (12.257, 57.03, 3.31), "Segin": (1.907, 63.67, 3.38),
"Meissa": (5.585, 9.93, 3.39),
}
FIGURES = [
["Betelgeuse", "Bellatrix", "Mintaka", "Rigel"],
["Mintaka", "Alnilam", "Alnitak", "Saiph"],
["Betelgeuse", "Alnitak"], ["Bellatrix", "Meissa", "Betelgeuse"],
["Dubhe", "Merak", "Phecda", "Megrez", "Dubhe"],
["Megrez", "Alioth", "Mizar", "Alkaid"],
["Segin", "Ruchbah", "Gamma Cas", "Schedar", "Caph"],
["Acrux", "Gacrux"], ["Mimosa", "Delta Cru"],
["Castor", "Pollux"], ["Vega", "Deneb", "Altair", "Vega"],
]
TEX_W, TEX_H = 1440, 720
NGP_RA, NGP_DEC = np.radians(192.85948), np.radians(27.12825)
GC_RA, GC_DEC = np.radians(266.41684), np.radians(-29.00781)
def sky_pixel(ra_h, dec_deg):
"""Catalogue coordinates → (column, row) in the equirectangular map."""
return (ra_h / 24.0 * (TEX_W - 1), (90.0 - dec_deg) / 180.0 * (TEX_H - 1))
def angsep(ra, dec, ra0, dec0):
return np.arccos(np.clip(
np.sin(dec) * np.sin(dec0) + np.cos(dec) * np.cos(dec0) * np.cos(ra - ra0),
-1.0, 1.0))
def _splat(img, x, y, radius, rgb, amp):
r = int(np.ceil(3 * radius))
xs = np.arange(int(round(x)) - r, int(round(x)) + r + 1)
ys = np.arange(int(round(y)) - r, int(round(y)) + r + 1)
ys = ys[(ys >= 0) & (ys < img.shape[0])]
if len(ys) == 0:
return
g = np.exp(-((ys[:, None] - y) ** 2 + (xs[None, :] - x) ** 2)
/ (2 * radius ** 2)) * amp
img[np.ix_(ys, xs % img.shape[1])] += g[..., None] * np.asarray(rgb)
def _stroke(img, p0, p1, rgb, amp):
(x0, y0), (x1, y1) = p0, p1
if abs(x1 - x0) > img.shape[1] / 2: # take the short way in RA
x1 += img.shape[1] if x0 > x1 else -img.shape[1]
n = int(max(abs(x1 - x0), abs(y1 - y0))) * 2 + 2
xs = np.linspace(x0, x1, n).astype(int) % img.shape[1]
ys = np.clip(np.linspace(y0, y1, n), 0, img.shape[0] - 1).astype(int)
img[ys, xs] = np.maximum(img[ys, xs], np.asarray(rgb) * amp)
# ── Build the sky image ──────────────────────────────────────────────────────
rng = np.random.default_rng(20260730)
ra_grid = np.linspace(0, 2 * np.pi, TEX_W)[None, :]
dec_grid = np.radians(np.linspace(90, -90, TEX_H))[:, None]
sky = np.zeros((TEX_H, TEX_W, 3), np.float32)
sky[:] = np.array([0.020, 0.028, 0.062], np.float32)
sky += (np.abs(np.sin(dec_grid)) ** 3 * 0.02)[..., None]
gal_lat = np.arcsin(np.sin(dec_grid) * np.sin(NGP_DEC)
+ np.cos(dec_grid) * np.cos(NGP_DEC)
* np.cos(ra_grid - NGP_RA))
band = np.exp(-(gal_lat / np.radians(11.0)) ** 2)
band = band * (0.45 + 0.85 * np.exp(
-(angsep(ra_grid, dec_grid, GC_RA, GC_DEC) / np.radians(55.0)) ** 2))
dust, amp = np.zeros((TEX_H, TEX_W), np.float32), 1.0
for k in range(1, 6):
dust += amp * (np.sin(3 * k * ra_grid + rng.uniform(0, 2 * np.pi))
* np.sin(2.5 * k * (dec_grid + np.pi / 2)
+ rng.uniform(0, 2 * np.pi)))
amp *= 0.55
band = band * np.clip(0.72 + 0.28 * dust, 0.25, 1.0)
sky += band[..., None] * np.array([0.46, 0.44, 0.55], np.float32)
n_faint = 24_000
f_dec = np.degrees(np.arcsin(rng.uniform(-1, 1, n_faint)))
fx = rng.uniform(0, TEX_W, n_faint).astype(int) % TEX_W
fy = ((90.0 - f_dec) / 180.0 * (TEX_H - 1)).astype(int)
keep = rng.random(n_faint) < (0.22 + 0.78 * band[fy, fx])
fx, fy = fx[keep], fy[keep]
sky[fy, fx] += rng.uniform(0.25, 0.95, (len(fx), 1)).astype(np.float32)
for chain in FIGURES:
for a, b in zip(chain, chain[1:]):
_stroke(sky, sky_pixel(*STARS[a][:2]), sky_pixel(*STARS[b][:2]),
(0.30, 0.58, 0.85), 1.0)
for name, (ra_h, dec_deg, mag) in STARS.items():
x, y = sky_pixel(ra_h, dec_deg)
flux = 10 ** (-0.4 * mag)
_splat(sky, x, y, 1.3 + 2.0 * flux ** 0.35, (1.0, 0.97, 0.90),
0.85 + 2.2 * flux ** 0.30)
sky_u8 = (np.clip(sky, 0, 1) * 255).astype(np.uint8)
# ── Spherical geometry ───────────────────────────────────────────────────────
NDEC, NRA = 49, 97
_ra = np.linspace(0, 2 * np.pi, NRA)
_dec = np.linspace(np.pi / 2, -np.pi / 2, NDEC)
RA, DEC = np.meshgrid(_ra, _dec)
X, Y, Z = np.cos(DEC) * np.cos(RA), np.cos(DEC) * np.sin(RA), np.sin(DEC)
def camera_to_zenith(azimuth, elevation):
"""Globe camera angles → the sky point facing you, as (RA hours, Dec°).
The renderer's turntable faces the unit vector with ``el = asin(vz)`` and
``az = atan2(vx, -vy)``. Feeding that through this example's
``x = cos δ cos α`` parameterisation collapses to ``α = azimuth − 90°``.
"""
return ((azimuth - 90.0) % 360.0) / 15.0, elevation
def zenith_to_camera(lst_hours, latitude):
"""The inverse: put (LST, latitude) overhead."""
return (lst_hours * 15.0 + 90.0) % 360.0, latitude
def horizon_segments(lst_hours, latitude, n=240):
"""Map-space polyline for the horizon: the circle 90° from the zenith.
Built as an orthonormal frame about the zenith rather than by solving the
altitude equation, so it stays well conditioned at the poles. Returned as
``(N, 2, 2)`` segments, split wherever the curve wraps in right ascension
so nothing is drawn straight across the map.
"""
a0, d0 = np.radians(lst_hours * 15.0), np.radians(latitude)
n_hat = np.array([np.cos(d0) * np.cos(a0), np.cos(d0) * np.sin(a0), np.sin(d0)])
# Any vector not parallel to the zenith gives a usable first basis vector.
tmp = np.array([0.0, 0.0, 1.0]) if abs(n_hat[2]) < 0.9 else np.array([1.0, 0.0, 0.0])
e1 = np.cross(n_hat, tmp)
e1 /= np.linalg.norm(e1)
e2 = np.cross(n_hat, e1)
t = np.linspace(0, 2 * np.pi, n)
pts = np.cos(t)[:, None] * e1 + np.sin(t)[:, None] * e2
ra_h = (np.degrees(np.arctan2(pts[:, 1], pts[:, 0])) % 360.0) / 15.0
dec_d = np.degrees(np.arcsin(np.clip(pts[:, 2], -1, 1)))
px = np.column_stack(sky_pixel(ra_h, dec_d))
segs = np.stack([px[:-1], px[1:]], axis=1) # (N-1, 2, 2)
keep = np.abs(segs[:, 0, 0] - segs[:, 1, 0]) < TEX_W / 2 # drop wrap jumps
return segs[keep]
def visible_stars(lst_hours, latitude):
"""Map pixels of the catalogue stars currently above the horizon."""
a0, d0 = np.radians(lst_hours * 15.0), np.radians(latitude)
out = []
for ra_h, dec_deg, _ in STARS.values():
alt = np.pi / 2 - angsep(np.radians(ra_h * 15.0), np.radians(dec_deg), a0, d0)
if alt > 0:
out.append(sky_pixel(ra_h, dec_deg))
return np.asarray(out) if out else np.zeros((0, 2))
# ── Figure ───────────────────────────────────────────────────────────────────
LST0, LAT0 = 5.5, 35.0 # Orion overhead, mid-northern latitude
fig = apl.Figure(figsize=(1080, 620))
gs = apl.GridSpec(3, 2, width_ratios=[1, 1.35], height_ratios=[6, 1, 1])
az0, el0 = zenith_to_camera(LST0, LAT0)
globe = fig.add_subplot(gs[0, 0]).plot_surface(
X, Y, Z, bounds=((-1, 1),) * 3, azimuth=az0, elevation=el0)
globe.set_axis_off()
globe.set_texture(sky_u8, cull_backfaces=True)
globe.set_title("celestial sphere — drag to move your zenith")
sky_map = fig.add_subplot(gs[0, 1]).imshow(sky_u8)
sky_map.set_title("sky map — inside the curve is above the horizon")
clock = fig.add_subplot(gs[1, :]).plot(
np.zeros(2), axes=[np.array([0.0, 24.0])], color="#37474f", linewidth=1)
clock.set_xlabel("local sidereal time (h)")
lat_ax = fig.add_subplot(gs[2, :]).plot(
np.zeros(2), axes=[np.array([-90.0, 90.0])], color="#37474f", linewidth=1)
lat_ax.set_xlabel("observer latitude ($\\degree$)")
# ── Overlays, redrawn whenever the zenith moves ──────────────────────────────
horizon = sky_map.add_lines(horizon_segments(LST0, LAT0), name="horizon",
edgecolors="#4dd0e1", linewidths=2.0)
# size_units="px" keeps the markers a fixed screen size as the map is zoomed.
zenith = sky_map.add_points(np.array([sky_pixel(LST0, LAT0)]), name="zenith",
sizes=9, color="#ffd54f", facecolors="#ffd54f",
alpha=1.0, size_units="px")
up_now = sky_map.add_points(visible_stars(LST0, LAT0), name="up",
sizes=4, color="#4dd0e1", facecolors="#4dd0e1",
alpha=0.55, size_units="px")
lst_line = clock.add_vline_widget(LST0, color="#ffd54f")
lat_line = lat_ax.add_vline_widget(LAT0, color="#ffd54f")
_syncing = False # guards the two-way coupling against feedback
def refresh(lst_hours, latitude, *, move_globe, move_sliders):
"""Push a new zenith everywhere. One writer, three followers."""
global _syncing
if _syncing:
return
_syncing = True
try:
with fig.batch():
horizon.set(segments=horizon_segments(lst_hours, latitude))
zenith.set(offsets=np.array([sky_pixel(lst_hours, latitude)]))
up_now.set(offsets=visible_stars(lst_hours, latitude))
sky_map.set_title(
f"sky map — LST {lst_hours:04.1f} h, latitude {latitude:+.0f}°")
if move_globe:
az, el = zenith_to_camera(lst_hours, latitude)
globe.set_view(azimuth=az, elevation=el)
if move_sliders:
lst_line.set(x=lst_hours)
lat_line.set(x=latitude)
finally:
_syncing = False
@globe.add_event_handler("pointer_move")
def _on_orbit(event):
"""Orbiting the globe IS choosing a time and a place."""
if event.azimuth is None:
return
lst_hours, latitude = camera_to_zenith(event.azimuth, event.elevation)
refresh(lst_hours, latitude, move_globe=False, move_sliders=True)
@lst_line.add_event_handler("pointer_move")
def _on_clock(event):
refresh(lst_line.x % 24.0, lat_line.x, move_globe=True, move_sliders=False)
@lat_line.add_event_handler("pointer_move")
def _on_latitude(event):
refresh(lst_line.x % 24.0, np.clip(lat_line.x, -89, 89),
move_globe=True, move_sliders=False)
refresh(LST0, LAT0, move_globe=False, move_sliders=False)
fig
Reading it#
The cyan circle is the horizon and the gold dot is the zenith. Stars inside the circle are above the ground and are marked in cyan; the rest have set.
Two things fall out of the geometry that are worth watching for:
Push the latitude to \(+90°\) and the horizon collapses onto the celestial equator — at the north pole the sky never rises or sets, it just wheels around Polaris.
Run the clock through 24 h at a mid latitude and the horizon sweeps the whole sky except two caps around the poles. Those caps are the circumpolar stars, the ones that never set from where you are standing.
refresh(2.0, 89.0, move_globe=True, move_sliders=True)
fig
Total running time of the script: (0 minutes 2.497 seconds)