Files
resolved/scripts/fetch-stars.py
T
2026-10-16 14:43:50 -04:00

400 lines
19 KiB
Python
Executable File

#!/usr/bin/env python3
"""Downloads the Hipparcos, Tycho-2 and Gaia DR3 star catalogs and writes the renderer's
star file.
Sources:
- Hipparcos main catalog, ESA 1997 (CDS I/239, hip_main.dat): 118,218 entries,
complete to about V 7.3, with Johnson V and B-V. Positions are ICRS at epoch
J1991.25. The 2007 new reduction (I/311) improves the astrometry by about a
milliarcsecond, far below any pixel here, but carries no V, so the original is used.
- Tycho-2, Hog et al. 2000 (CDS I/259, tyc2.dat.00-19): 2.54 million stars to about
VT 11.5 (99 % complete to V 11.0), with Tycho BT and VT magnitudes. Mean positions
are ICRS at epoch J2000.
- Gaia DR3, Gaia Collaboration 2023 (gaiadr3.gaia_source_lite, from the ESA Gaia
archive's TAP service, https://gea.esac.esa.int/archive/): every star brighter than
G 15 by default (--gaia), 36.9 million, with its G magnitude and BP - RP color.
Positions are ICRS at epoch J2016.0. Gaia's photometry saturates for the brightest
stars, which Hipparcos and Tycho-2 cover.
Every star is taken from one catalog: Hipparcos first, then Tycho-2 for the stars it
does not have (by Hipparcos number), then Gaia for the stars neither has (by Gaia's own
cross-matches to them, gaiadr3.hipparcos2_best_neighbour and
gaiadr3.tycho2tdsc_merge_best_neighbour). Magnitudes are converted to Johnson V and
B - V:
- Tycho: V = VT - 0.090 (BT - VT), B - V = 0.850 (BT - VT) (ESA 1997, vol. 1, 1.3,
appendix 4). A star with no BT keeps V = VT and no color.
- Gaia: V from G and BP - RP by Riello et al. 2021 (A&A 649, A3, table C.2),
G - V = -0.02704 + 0.01424 c - 0.2156 c^2 + 0.01426 c^3 for c = BP - RP in -0.5 to
2.75, 0.03 mag scatter. B - V from BP - RP by a cubic fitted here to the Hipparcos
stars Gaia matched (printed when the script runs). A star with no BP - RP gets
c = 0.82, the median of stars near G 15, and no B - V.
Positions are carried to J2000 by their proper motion. The few Tycho-2 stars with no
mean position (pflag X) use their observed position, at about epoch 1991.5, with no
proper motion. Parallax is not kept: it moves every star by under 1 arcsecond.
Hipparcos and Tycho-2 come from the CDS archive (https://cdsarc.cds.unistra.fr) and are
checked by sha256. Gaia comes from about 2,000 queries of the archive, each a range of
source_id (a patch of sky), 4 at a time, split in half when a reply reaches the
service's 50,000-row limit; a query's results are kept, so a run can be resumed, but
they are not pinned by checksum. At G 15 that is about 2.3 GB of text and takes about
an hour. Everything goes to data/stars/source/; the output is data/stars/stars.bin,
about 600 MB at G 15, both gitignored. Needs numpy and pandas.
The output's format, little-endian:
8 bytes b"RSTARS02"
u32 number of stars
u32 0
then per star, brightest first, 16 bytes:
u32 right ascension, ICRS, epoch J2000, in units of 2 pi / 2^32
i32 declination, ICRS, epoch J2000, in units of pi / 2^32
i16 proper motion in right ascension times cos(declination), mas/yr
i16 proper motion in declination, mas/yr
u16 Johnson V magnitude plus 5, in thousandths
i16 Johnson B - V in thousandths, or -32768 if unknown
Usage: scripts/fetch-stars.py [--gaia G] [--faint V]
--gaia G take Gaia stars brighter than G (default 15; 0 for none)
--faint V keep stars to V (default: all those the catalogs give)
"""
import argparse
import concurrent.futures
import gzip
import hashlib
import io
import math
import os
import sys
import time
import urllib.parse
import urllib.request
import numpy as np
import pandas as pd
CDS = "https://cdsarc.cds.unistra.fr/ftp/I"
GAIA_TAP = "https://gea.esac.esa.int/tap-server/tap/sync"
# The service's limit on rows per synchronous query.
GAIA_ROWS = 50_000
# Gaia source_id carries the HEALPix level-12 cell in its high bits.
CELLS = 12 * 4**12
CELL = 2**35
FILES = {
"hip_main.dat": ("239/hip_main.dat", "58ceabb104d647160d9437ce6e513a02a036bb4ad9f8879a5a22fd52943616e0"),
"tyc2.dat.00.gz": ("259/tyc2.dat.00.gz", "1b224570fc6eb151984ce106fdf797728649d1ea77aaf3effd2d3444cfac6df6"),
"tyc2.dat.01.gz": ("259/tyc2.dat.01.gz", "36221f4a5cdd9c5009d64299d2e6409f8d4cdf924321768780165bd3e2a2a99a"),
"tyc2.dat.02.gz": ("259/tyc2.dat.02.gz", "e2ea4eeb7d204dde70ac54d172b1c429540c0d23d62ec0d7500561adc44cc57b"),
"tyc2.dat.03.gz": ("259/tyc2.dat.03.gz", "29e0cb57cba7e1651455efcb528276b0c28e1fe83f47f8ebe963e48d8b180afb"),
"tyc2.dat.04.gz": ("259/tyc2.dat.04.gz", "d04868329b470aa320f0fe5e9c525be619a7a228b364433e5f114f2178946d9b"),
"tyc2.dat.05.gz": ("259/tyc2.dat.05.gz", "51f08772197dcf0cd26ff098c866036340b797f71ee04be22712517c10c73c4e"),
"tyc2.dat.06.gz": ("259/tyc2.dat.06.gz", "f55503a66abce0e11d1a4ccbcde1f4c97a6e3e769de85e005f57ba963086ab03"),
"tyc2.dat.07.gz": ("259/tyc2.dat.07.gz", "bc15b5b1f6308477360fbb93f6a98d9029601cab68973202be0f8764c127f524"),
"tyc2.dat.08.gz": ("259/tyc2.dat.08.gz", "315ac34c678bc2e0f568b43b6b2d7c3d7c9cf9c089b51115d18b872864e50426"),
"tyc2.dat.09.gz": ("259/tyc2.dat.09.gz", "d66c671fa29aad10bc5fd697f68c918bb774c1ffc3c7db9c20234895268390d5"),
"tyc2.dat.10.gz": ("259/tyc2.dat.10.gz", "d86529df819ebdf3f4b8892510973ee10e1db06930814e014dd467193e740d52"),
"tyc2.dat.11.gz": ("259/tyc2.dat.11.gz", "fc0c3203b3d2787da54c43d9963e363c8a3926576c383b101a15350a7ca23e9c"),
"tyc2.dat.12.gz": ("259/tyc2.dat.12.gz", "68fd9d7e7353d52dea0de043cc62b45c04d234104eb05c063dc23ea7b6759576"),
"tyc2.dat.13.gz": ("259/tyc2.dat.13.gz", "a08ecc092e6d134742c0f594c56c5c2e02b8fd5e76cffb1779877f2fcefac3fc"),
"tyc2.dat.14.gz": ("259/tyc2.dat.14.gz", "2846e8cc489795a4d8160d2d5b64ff0e9a00598adba35a52a4f8ccf4c5e38b3e"),
"tyc2.dat.15.gz": ("259/tyc2.dat.15.gz", "74784a3656382090a70d778820335f4a781509dbe3fcb9e92d18cf96b2c46c72"),
"tyc2.dat.16.gz": ("259/tyc2.dat.16.gz", "320560e3d551cc40fa1f54ef7133709bf4bd45efa6688f808be7c09be5544a4d"),
"tyc2.dat.17.gz": ("259/tyc2.dat.17.gz", "a8154e940aa0a0e8f31d91b4bd9fd56e6849da42cddb13c7dda772439b06991a"),
"tyc2.dat.18.gz": ("259/tyc2.dat.18.gz", "5e951a1a51df7956f205b83cbfa5501d357c843640a253c1e6d3917bebe7d928"),
"tyc2.dat.19.gz": ("259/tyc2.dat.19.gz", "f59605a38116f517a31a7dbdee3469c077658f2f40b8afe5da2aeb832eaee3dd"),
}
# A star's fields as arrays: right ascension and declination (radians, J2000), proper
# motion (mas/yr), V and B - V (NaN if unknown).
FIELDS = ("ra", "dec", "pm_ra", "pm_dec", "v", "bv")
def fetch(source):
os.makedirs(source, exist_ok=True)
for name, (path, digest) in FILES.items():
target = os.path.join(source, name)
if not os.path.exists(target):
print(f"fetching {name}")
urllib.request.urlretrieve(f"{CDS}/{path}", target + ".part")
os.replace(target + ".part", target)
h = hashlib.sha256()
with open(target, "rb") as f:
for block in iter(lambda: f.read(1 << 20), b""):
h.update(block)
if h.hexdigest() != digest:
sys.exit(f"{target}: sha256 {h.hexdigest()}, expected {digest}")
def number(line, first, last):
"""The value in 1-based columns first to last, or None if blank."""
text = line[first - 1:last].strip()
return float(text) if text else None
def propagate(ra, dec, pm_ra, pm_dec, years):
"""Moves positions (radians) by their proper motion (mas/yr, right ascension times
cos declination) over `years`, along the sphere."""
scale = np.radians(1 / 3.6e6) * years
ca, sa, cd, sd = np.cos(ra), np.sin(ra), np.cos(dec), np.sin(dec)
x = cd * ca + scale * (pm_ra * -sa + pm_dec * -sd * ca)
y = cd * sa + scale * (pm_ra * ca + pm_dec * -sd * sa)
z = sd + scale * (pm_dec * cd)
norm = np.sqrt(x * x + y * y + z * z)
return np.mod(np.arctan2(y, x), 2 * np.pi), np.arcsin(z / norm)
def table(rows):
columns = list(zip(*rows)) if rows else [[]] * 6
return {k: np.array(c, dtype=float) for k, c in zip(FIELDS, columns)}
def hipparcos(path):
rows, numbers = [], set()
with open(path, encoding="ascii") as f:
for line in f:
v, ra, dec = number(line, 42, 46), number(line, 52, 63), number(line, 65, 76)
if v is None or ra is None or dec is None:
continue
bv = number(line, 246, 251)
rows.append((math.radians(ra), math.radians(dec), number(line, 88, 95) or 0.0,
number(line, 97, 104) or 0.0, v, math.nan if bv is None else bv))
numbers.add(int(line[8:14]))
stars = table(rows)
stars["ra"], stars["dec"] = propagate(
stars["ra"], stars["dec"], stars["pm_ra"], stars["pm_dec"], 2000.0 - 1991.25)
return stars, numbers
def hipparcos_colors(path):
"""B - V by Hipparcos number."""
colors = {}
with open(path, encoding="ascii") as f:
for line in f:
bv = number(line, 246, 251)
if bv is not None:
colors[int(line[8:14])] = bv
return colors
def tycho2(paths, hip_numbers):
rows, skipped = [], {"hipparcos": 0, "no VT": 0}
for path in paths:
with gzip.open(path, "rt", encoding="ascii") as f:
for line in f:
hip = line[142:148].strip()
if hip and int(hip) in hip_numbers:
skipped["hipparcos"] += 1
continue
bt, vt = number(line, 111, 116), number(line, 124, 129)
if vt is None:
skipped["no VT"] += 1
continue
if line[13] == "X":
ra, dec = number(line, 153, 164), number(line, 166, 177)
pm_ra = pm_dec = 0.0
else:
ra, dec = number(line, 16, 27), number(line, 29, 40)
pm_ra, pm_dec = number(line, 42, 48) or 0.0, number(line, 50, 56) or 0.0
if bt is None:
v, bv = vt, math.nan
else:
v, bv = vt - 0.090 * (bt - vt), 0.850 * (bt - vt)
rows.append((math.radians(ra), math.radians(dec), pm_ra, pm_dec, v, bv))
return table(rows), skipped
def gaia_query(query):
"""The rows of an ADQL query as a DataFrame, retrying a few times."""
data = urllib.parse.urlencode({"REQUEST": "doQuery", "LANG": "ADQL", "FORMAT": "csv",
"QUERY": query}).encode()
for attempt in range(6):
try:
with urllib.request.urlopen(GAIA_TAP, data=data, timeout=600) as reply:
return pd.read_csv(io.BytesIO(reply.read()))
except Exception as e: # network or service errors: wait and retry
if attempt == 5:
raise
print(f" retrying after {e}")
time.sleep(10 * (attempt + 1))
def gaia_range(directory, faint, lo, hi):
"""Gaia stars brighter than G `faint` with source_id in [lo, hi), from the cache in
`directory` or the archive, split until each query is within the row limit."""
path = os.path.join(directory, f"{lo}-{hi}.csv")
split = os.path.join(directory, f"{lo}-{hi}.split")
mid = (lo + hi) // 2
if os.path.exists(path):
return [path]
if os.path.exists(split):
return gaia_range(directory, faint, lo, mid) + gaia_range(directory, faint, mid, hi)
query = (
"SELECT g.ra, g.dec, g.pmra, g.pmdec, g.phot_g_mean_mag, g.bp_rp, "
"t.source_id AS tycho, h.source_id AS hip "
"FROM gaiadr3.gaia_source_lite AS g "
"LEFT OUTER JOIN gaiadr3.tycho2tdsc_merge_best_neighbour AS t ON g.source_id = t.source_id "
"LEFT OUTER JOIN gaiadr3.hipparcos2_best_neighbour AS h ON g.source_id = h.source_id "
f"WHERE g.phot_g_mean_mag < {faint} AND g.source_id >= {lo} AND g.source_id < {hi}"
)
rows = gaia_query(query)
if len(rows) >= GAIA_ROWS:
# Too many for one reply: split, and remember it, so a rerun does not ask again.
open(split, "w").close()
return gaia_range(directory, faint, lo, mid) + gaia_range(directory, faint, mid, hi)
rows.to_csv(path + ".part", index=False)
os.replace(path + ".part", path)
return [path]
def gaia(source, faint):
directory = os.path.join(source, f"gaia-g{faint:g}")
os.makedirs(directory, exist_ok=True)
chunks = 2048
edges = [CELLS * i // chunks * CELL for i in range(chunks + 1)]
paths = []
with concurrent.futures.ThreadPoolExecutor(4) as pool:
jobs = [pool.submit(gaia_range, directory, faint, lo, hi) for lo, hi in zip(edges, edges[1:])]
for done, job in enumerate(concurrent.futures.as_completed(jobs), 1):
paths += job.result()
if done % 100 == 0:
print(f" Gaia: {done} of {chunks} ranges")
return pd.concat((pd.read_csv(p) for p in sorted(paths)), ignore_index=True)
def gaia_color_fit(source, hip_main):
"""A cubic giving B - V from BP - RP, fitted to the Hipparcos stars Gaia matched, and
the scatter of the fit."""
path = os.path.join(source, "gaia-hipparcos.csv")
if not os.path.exists(path):
parts = []
edges = [CELLS * i // 8 * CELL for i in range(9)]
for lo, hi in zip(edges, edges[1:]):
parts.append(gaia_query(
"SELECT h.original_ext_source_id AS hip, g.bp_rp "
"FROM gaiadr3.hipparcos2_best_neighbour AS h "
"JOIN gaiadr3.gaia_source_lite AS g ON h.source_id = g.source_id "
f"WHERE g.source_id >= {lo} AND g.source_id < {hi}"
))
pd.concat(parts).to_csv(path, index=False)
matched = pd.read_csv(path).dropna()
colors = hipparcos_colors(hip_main)
c = matched["bp_rp"].to_numpy()
bv = np.array([colors.get(int(h), np.nan) for h in matched["hip"]])
ok = np.isfinite(bv) & (c > -0.5) & (c < 3.5)
coefficients = np.polyfit(c[ok], bv[ok], 3)
residual = bv[ok] - np.polyval(coefficients, c[ok])
return coefficients, np.median(np.abs(residual)) * 1.4826, ok.sum()
def gaia_stars(rows, coefficients):
"""Gaia rows that neither Hipparcos nor Tycho-2 has, as V and B - V at J2000."""
rows = rows[rows["tycho"].isna() & rows["hip"].isna()]
c = rows["bp_rp"].to_numpy()
known = np.isfinite(c)
c = np.where(known, c, 0.82)
cc = np.clip(c, -0.5, 2.75)
g_minus_v = -0.02704 + 0.01424 * cc - 0.2156 * cc**2 + 0.01426 * cc**3
pm_ra = np.nan_to_num(rows["pmra"].to_numpy())
pm_dec = np.nan_to_num(rows["pmdec"].to_numpy())
ra, dec = propagate(np.radians(rows["ra"].to_numpy()), np.radians(rows["dec"].to_numpy()),
pm_ra, pm_dec, 2000.0 - 2016.0)
return {
"ra": ra, "dec": dec, "pm_ra": pm_ra, "pm_dec": pm_dec,
"v": rows["phot_g_mean_mag"].to_numpy() - g_minus_v,
"bv": np.where(known, np.polyval(coefficients, np.clip(c, -0.5, 3.5)), np.nan),
}
def unit(stars):
cd = np.cos(stars["dec"])
return np.stack([cd * np.cos(stars["ra"]), cd * np.sin(stars["ra"]), np.sin(stars["dec"])], 1)
def drop_near(stars, others, radius):
"""`stars` without those within `radius` radians of any of `others`, and how many
were dropped: cells of the radius on a side, searched with their neighbors."""
key = lambda c: ((c[:, 0] + (1 << 20)) << 42) | ((c[:, 1] + (1 << 20)) << 21) | (c[:, 2] + (1 << 20))
points = unit(stars)
keys = key(np.floor(points / radius).astype(np.int64))
order = np.argsort(keys)
# A key no cell has ends the list, so a walk past the last star stops there.
sorted_keys = np.append(keys[order], np.iinfo(np.int64).max)
near = unit(others)
cells = np.floor(near / radius).astype(np.int64)
drop = np.zeros(len(points), dtype=bool)
for offset in np.ndindex(3, 3, 3):
k = key(cells + np.array(offset) - 1)
i = np.searchsorted(sorted_keys, k)
# Every star in the cell: walk forward while the key repeats.
while True:
hit = sorted_keys[i] == k
if not hit.any():
break
j = order[i[hit]]
close = np.linalg.norm(points[j] - near[hit], axis=1) < radius
drop[j[close]] = True
i = i + hit
return {f: stars[f][~drop] for f in FIELDS}, int(drop.sum())
def write(path, stars, faint):
keep = stars["v"] <= faint
order = np.argsort(stars["v"][keep], kind="stable")
s = {k: stars[k][keep][order] for k in FIELDS}
record = np.dtype([("ra", "<u4"), ("dec", "<i4"), ("pm_ra", "<i2"), ("pm_dec", "<i2"),
("v", "<u2"), ("bv", "<i2")])
out = np.zeros(len(s["v"]), dtype=record)
out["ra"] = np.round(np.mod(s["ra"], 2 * np.pi) / (2 * np.pi) * 2**32).astype(np.uint64) % 2**32
out["dec"] = np.clip(np.round(s["dec"] / np.pi * 2**32), -(2**31), 2**31 - 1)
out["pm_ra"] = np.clip(np.round(s["pm_ra"]), -32767, 32767)
out["pm_dec"] = np.clip(np.round(s["pm_dec"]), -32767, 32767)
out["v"] = np.clip(np.round((s["v"] + 5.0) * 1000), 0, 65535)
out["bv"] = np.where(np.isfinite(s["bv"]), np.clip(np.round(s["bv"] * 1000), -32767, 32767), -32768)
with open(path + ".part", "wb") as f:
f.write(b"RSTARS02" + np.array([len(out), 0], dtype="<u4").tobytes())
f.write(out.tobytes())
os.replace(path + ".part", path)
return len(out), s["v"][-1] if len(out) else math.nan
def main():
parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
parser.add_argument("--gaia", type=float, default=15.0, help="Gaia stars brighter than G")
parser.add_argument("--faint", type=float, default=math.inf, help="faintest V kept")
args = parser.parse_args()
root = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "data", "stars")
source = os.path.join(root, "source")
fetch(source)
hip_main = os.path.join(source, "hip_main.dat")
hip, numbers = hipparcos(hip_main)
tyc, skipped = tycho2(
sorted(os.path.join(source, n) for n in FILES if n.startswith("tyc2")), numbers
)
print(f"Hipparcos: {len(hip['v'])} stars; Tycho-2: {len(tyc['v'])} more "
f"({skipped['hipparcos']} in Hipparcos, {skipped['no VT']} without VT)")
parts = [hip, tyc]
if args.gaia > 0:
coefficients, scatter, used = gaia_color_fit(source, hip_main)
print(f"B - V from BP - RP: cubic {np.round(coefficients, 5).tolist()} "
f"(highest power first), {scatter:.3f} mag scatter over {used} Hipparcos stars")
rows = gaia(source, args.gaia)
extra = gaia_stars(rows, coefficients)
print(f"Gaia DR3 to G {args.gaia:g}: {len(rows)} stars, {len(extra['v'])} not in "
f"Hipparcos or Tycho-2 by Gaia's cross-match")
# Gaia's cross-match misses about 1 % of Tycho-2: drop Gaia stars within 1
# arcsecond of one, which Tycho's and Hipparcos's photometry, unable to resolve
# pairs that close, already includes.
others = {k: np.concatenate([hip[k], tyc[k]]) for k in FIELDS}
extra, dropped = drop_near(extra, others, np.radians(1 / 3600))
print(f" {dropped} more within 1 arcsecond of a Hipparcos or Tycho-2 star")
parts.append(extra)
stars = {k: np.concatenate([p[k] for p in parts]) for k in FIELDS}
out = os.path.join(root, "stars.bin")
count, faintest = write(out, stars, args.faint)
print(f"{count} stars to V {faintest:.2f} in {out}")
if __name__ == "__main__":
main()