125 lines
3.8 KiB
Python
125 lines
3.8 KiB
Python
"""Cross-check the ApsisRS geodetic to Earth-fixed conversion.
|
||
|
||
Reads the grid written by ``cargo run -p apsis-rs --features physics
|
||
--example dump_geodesy`` and recomputes each position with astropy, then
|
||
reports the largest difference.
|
||
|
||
The script needs no network access. It reads the input from a file path
|
||
given on the command line, defaulting to the standard location, and writes a
|
||
report to standard output.
|
||
|
||
Usage:
|
||
|
||
pixi run python validation/python/coord/check_geodetic.py [path]
|
||
|
||
Exit status is zero when every position agrees within the tolerance and one
|
||
otherwise.
|
||
|
||
交叉核对 ApsisRS 的大地坐标到地固坐标转换。
|
||
|
||
读取由 ``cargo run -p apsis-rs --features physics --example dump_geodesy``
|
||
写出的网格，并用 astropy 逐个重算位置，然后报告最大差异。
|
||
|
||
本脚本不需要网络访问。它从命令行给出的文件路径读取输入，缺省为标准位置，
|
||
并把报告写到标准输出。
|
||
|
||
用法：
|
||
|
||
pixi run python validation/python/coord/check_geodetic.py [path]
|
||
|
||
每个位置都在容差内一致时退出状态为零，否则为一。
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import csv
|
||
import sys
|
||
from pathlib import Path
|
||
|
||
import astropy.units as u
|
||
from astropy.coordinates import EarthLocation
|
||
|
||
# Position tolerance in meters. Astropy evaluates the same closed form, so
|
||
# the only difference should be floating point rounding.
|
||
# 位置容差，单位为米。Astropy 求的是同一个闭式，故唯一差异应是浮点舍入。
|
||
TOLERANCE_M = 1e-6
|
||
|
||
DEFAULT_INPUT = Path("validation/data/expected/rust_geodetic.csv")
|
||
|
||
|
||
def rust_positions(path: Path) -> list[tuple[float, float, float, tuple[float, float, float]]]:
|
||
"""Reads the Rust output rows.
|
||
|
||
读取 Rust 输出各行。
|
||
"""
|
||
rows = []
|
||
with path.open(newline="") as handle:
|
||
reader = csv.DictReader(handle)
|
||
for row in reader:
|
||
rows.append(
|
||
(
|
||
float(row["latitude_degrees"]),
|
||
float(row["longitude_degrees"]),
|
||
float(row["height_m"]),
|
||
(float(row["x_m"]), float(row["y_m"]), float(row["z_m"])),
|
||
)
|
||
)
|
||
return rows
|
||
|
||
|
||
def astropy_position(latitude_deg: float, longitude_deg: float, height_m: float):
|
||
"""Recomputes a position with astropy.
|
||
|
||
用 astropy 重算一个位置。
|
||
"""
|
||
location = EarthLocation.from_geodetic(
|
||
lon=longitude_deg * u.deg,
|
||
lat=latitude_deg * u.deg,
|
||
height=height_m * u.m,
|
||
)
|
||
return (
|
||
location.x.to_value(u.m),
|
||
location.y.to_value(u.m),
|
||
location.z.to_value(u.m),
|
||
)
|
||
|
||
|
||
def main(argv: list[str]) -> int:
|
||
path = Path(argv[1]) if len(argv) > 1 else DEFAULT_INPUT
|
||
if not path.exists():
|
||
print(f"missing input file: {path}", file=sys.stderr)
|
||
print(
|
||
"generate it with: cargo run -p apsis-rs --features physics "
|
||
"--example dump_geodesy > " + str(path),
|
||
file=sys.stderr,
|
||
)
|
||
return 2
|
||
|
||
rows = rust_positions(path)
|
||
if not rows:
|
||
print("input file has no data rows", file=sys.stderr)
|
||
return 2
|
||
|
||
worst = 0.0
|
||
worst_at = None
|
||
for latitude, longitude, height, rust in rows:
|
||
reference = astropy_position(latitude, longitude, height)
|
||
error = max(abs(a - b) for a, b in zip(rust, reference))
|
||
if error > worst:
|
||
worst = error
|
||
worst_at = (latitude, longitude, height)
|
||
|
||
print(f"rows compared: {len(rows)}")
|
||
print(f"worst error: {worst:.3e} m at lat/lon/height {worst_at}")
|
||
print(f"tolerance: {TOLERANCE_M:.3e} m")
|
||
|
||
if worst <= TOLERANCE_M:
|
||
print("result: pass")
|
||
return 0
|
||
print("result: FAIL")
|
||
return 1
|
||
|
||
|
||
if __name__ == "__main__":
|
||
raise SystemExit(main(sys.argv))
|