Compute get_lonlatalt altitude on WGS84 consistently - #247
Conversation
cb946b1 to
f274403
Compare
get_lonlatalt scaled the position by XKMPER (the WGS72 radius SGP4 propagates with) and the resulting altitude by A (the WGS84 radius), which inflated every altitude by A / XKMPER - 1: about 2.1 m for a low orbit and 13 m at geostationary altitude. The position is now scaled by A, so the geodetic conversion is entirely on WGS84 while SGP4 keeps its WGS72 constants.
f274403 to
15854be
Compare
mraspaud
left a comment
There was a problem hiding this comment.
Great catch! a 1mm precision is fantastic :)
Thanks a lot for the clean fix
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #247 +/- ##
=======================================
Coverage 92.17% 92.17%
=======================================
Files 19 19
Lines 4127 4154 +27
=======================================
+ Hits 3804 3829 +25
- Misses 323 325 +2
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Manny7717
left a comment
There was a problem hiding this comment.
Independent verification confirms the mixed-datum bug and the fix. Details:
Bug is real (executed). On base origin/main (0a154c3), get_position(normalize=True) divides the km position by XKMPER (6378.135, WGS72), but the geodetic reduction then scales altitude by A (6378.137, WGS84) — dimensionally inconsistent by A/XKMPER − 1 ≈ 3.1e-7. Old code returns 392.019534 km altitude for the ISS test epoch vs 392.017412 km under a consistent WGS84 reduction: a +2.12 m inflation (matches the predicted ~2.3 m at r/cosφ ≈ 7264 km). At GEO radius the same relative error is ~13.2 m (arithmetic check).
Fix is correct. Using normalize=False (km) and scaling by A makes the position, the flattening F, and the altitude scale all refer to the same WGS84 ellipsoid. I verified the new implementation is algebraically identical to the standard geodetic iteration (position in A-radii, c = 1/sqrt(1−e²sin²φ) = N/A): alt = A·(r_norm/cosφ − c) ≡ r_km/cosφ − N. Lon is bit-identical (ratio unchanged); lat shifts only ~4e-8 deg (mixed-datum artifact removal, sub-mm, still within the existing eps_deg tolerance). utc2local (the only in-repo consumer) uses lon only — unaffected.
Regression-proven. The PR's new test_sublonlat_altitude_is_geodetic_on_wgs84 (external 50-iteration geodetic reference, tight 1e-9 km tolerance) FAILS on base code with exactly the +2.12 m discrepancy and PASSES on head. Updated test_sublonlat expected value likewise fails on base, passes on head.
Full suite: 247 passed / 1 skipped on head (test_orbital.py 29/29). New code lines within ruff line-length 120; comments accurate. No competing/overlapping PR touching Orbital.get_lonlatalt (checked open pyorbital PRs; #87's recent merge restructured get_observer_look only).
Non-blocking nit only: the method docstring could state the altitude is the WGS84 ellipsoidal height (it currently links celestrak v02n03 without naming the datum), and vel remains unpacked-but-unused as before — neither worth holding the PR.
Orbital.get_lonlataltconverts the propagated position to geodetic coordinates on the WGS84 ellipsoid, but it takes the position fromget_position(normalize=True), which divides kilometers byXKMPER= 6378.135 km, the WGS72 equatorial radius SGP4 propagates with (Hoots and Roehrich, Spacetrack Report No. 3, 1980). The geodetic iteration then runs with the WGS84 flattening, and the altitude is scaled back withA= 6378.137 km, the WGS84 equatorial radius (NIMA TR8350.2). Dividing by one radius and multiplying by the other inflates every altitude byA / XKMPER - 1: about 2.1 m for the ISS and about 13 m at geostationary altitude. Longitude and latitude are unaffected, since they depend only on direction.The conversion now scales the kilometer position by
A, so the whole geodetic step is on WGS84. SGP4 itself still runs with the WGS72 constants the elements are generated with, as the propagator's reference implementation recommends (Vallado, Crawford, Hujsak and Kelso, "Revisiting Spacetrack Report #3", AIAA 2006-6753).The existing
test_sublonlatexpected altitude moves from 392.01953 km to 392.01741 km, and its altitude tolerance tightens to 1 mm. A new test converts the same positions to geodetic coordinates independently in kilometers on WGS84 and checks the altitude to a nanometre across an orbit, and checks that the previously inflated value no longer matches.