Add ENU/geodetic image transforms needed by the GMN contrail pipeline

[KeatonBuddy]

What this adds

The GMN contrail pipeline (GMNContrailPipeline → Janus) needs a set of ENU/geodetic
image transforms that currently exist only on the test-coordinate-transforms branch.
Rather than merge that whole branch, this ports forward only the surface the pipeline
actually imports, on top of prerelease.

RMS/Astrometry/ApplyAstrometry.py — xyHtToENUPP, enHtToXYPP, enuToXYPP,
geoToENUPP, geoToXYPP, ENHt0ToENHt1, xyToAltAzPP, plus the imageCenter and
_centreAltAz helpers they share.

RMS/Astrometry/CyFunctions.pyx — cyXYToAltAz, cyXYHttoENU_wgs84, cyGeoToENU,
cyGeoToXY_wgs84_iter, cyENUToXY_iter, cyENHtToXY_iter, cyENHt0ToENHt1 and their
three geodetic cdef helpers. Appended; nothing existing was modified.

RMS/Astrometry/GPUENHt.py — new module (optional downstream, guarded by
except ImportError).

RMS/Formats/Platepar.py — sets height_wgs84 from elev via EGM96 in
loadFromDict() and read(), as the source branch does. The ENU/geodetic transforms
need the ellipsoidal station height rather than the MSL elev.

RMS/GeoidHeightEGM96.py — taken from the source branch. This file had not been
modified on prerelease since 073d33a5 (2021-06-28), and the merge-base blob is
identical to prerelease's, so the branch's file is prerelease's plus dd8e1383
(2025-09-05). Same four functions on both sides; nothing dropped.

What this deliberately does not change

Every function prerelease already has keeps prerelease's implementation. In
particular rotationWrtHorizon and xyHt2Geo are byte-identical to prerelease,
even though the source branch reworked both — those rewrites are not included here.

One adaptation was unavoidable: the ported functions call
pyRefractionTrueToApparent / pyRefractionApparentToTrue instead of
refractionTrueToApparent / refractionApparentToTrue, because prerelease declares
the latter cdef, so they are not callable from Python. The py* forms are
prerelease's own cpdef wrappers around the same code.

Verification

• prerelease astrometry is bit-identical across xyToRaDecPP, raDecToXYPP,
  rotationWrtHorizon, rotationWrtHorizonToPosAngle, computeFOVSize, xyHt2Geo,
  fovEdgePolygon, getFOVSelectionRadius, correctVignetting
• No prerelease symbol is lost (AST scan over every .py); the diff is additive
• Cython rebuilds from a cleared ~/.pyxbld; SkyFit2 starts; core modules import
• Round-trips check out: XY→ENU→XY and geodetic→XY→geodetic agree to ~0.01 px, and
  ENHt0ToENHt1 preserves the line of sight when changing target height

Known gap — needs a decision

dd8e1383 changed the GeoidHeightEGM96 signatures (config → optional path) but
never updated the callers that pass a Config. Those still pass one where a path is
now expected, and raise when reached:

Location	Function
RMS/Astrometry/Conversions.py:516,528	addECEFVectortoLatLonEle
RMS/Astrometry/Conversions.py:577,578	getECEFVectorBetweenGeoPoints
RMS/Formats/CAL.py:119	writeCAL
RMS/Formats/FrameInterface.py:1697	FrameInterface.__init__

This is pre-existing on the source branch (the pipeline simply never reaches these
paths). They are left untouched here because they are prerelease code this branch
does not otherwise modify. The fix is small — pass
os.path.join(config.egm96_path, config.egm96_file_name), or let the argument default —
but whether to make GeoidHeightEGM96 accept a Config for backward compatibility, or
update these four call sites, is a maintainer call. Happy to do either.

Marked draft for that reason.

🤖 Generated with Claude Code

[dvida]

Thanks for the careful port and the detailed write-up. I built the PR head (8e6a3149) from scratch and compared the new transforms against the existing prerelease astrometry on two real GMN platepars (US05F1, radial7-odd, centre alt 69°; NZ000N, radial7-odd, centre alt 36°). The Cython kernels themselves are numerically sound, but the Python wrappers feed them the wrong rotation angle, so the ported transforms currently disagree with xyToRaDecPP by tens of arcminutes on an unrolled camera and by many degrees on a rolled one. There are also two runtime regressions for existing users. Details below; inline comments point at the exact lines.

Blocking

B1. Wrong rotation angle passed to the alt/az gnomonic (rot = rotationWrtHorizon(platepar) in all five wrappers).

The kernels need the position angle of the image +x axis in the local tangent plane. Prerelease rotationWrtHorizon returns atan2(dAlt, dAz) without cos(alt) on the azimuth step; it is a display quantity, not a tangent-plane angle. Reference: xyToRaDecPP → cyTrueRaDec2ApparentAltAz(…, refraction=False) (true alt/az, same epoch), 63-point grid over the full frame:

platepar	roll added to pos_angle_ref	max error, as shipped
US05F1 (alt 69°)	0	28′
US05F1	+45°	1104′ (18.4°)
NZ000N (alt 36°)	0	38′
NZ000N	+45°	285′

Fitting the rotation numerically brings the max error down to 0.09′–0.2′ across the FOV, so the Cython math is right and only the input is wrong. The best-fit rotation is exactly

rot = pos_angle_ref − q,   q = atan2(sin HA, tan(lat)·cos(dec) − sin(dec)·cos(HA))

with q the parallactic angle at the FOV centre (HA from JD2LST(pp.JD, pp.lon) and pp.RA_d; RA_d/dec_d are epoch-of-date in RMS, which is also why the cyraDec2AltAz centre call without precession is correct). This matches the fit to 0.0003°–0.0007° at rolls 0/45/90 on both platepars. I'd suggest a new analytic helper (e.g. rotationWrtHorizonTangentPlane(platepar)) used by the five wrappers, and leaving prerelease rotationWrtHorizon (26 call sites) alone. Note that the source branch's finite-difference tangent-plane version is not accurate enough either (it left 26′–224′ residuals in my tests), so "numerically identical to the source branch" is not a useful acceptance criterion here.

B2. Refraction applied twice to the FOV centre. Each wrapper does alt_centre = pyRefractionTrueToApparent(alt_centre) unconditionally, then the kernel applies refractionTrueToApparent(h0) again when refraction=True. When platepar.refraction is False the Python side still refracts once. Either way the centre carries one extra refraction: 0.47′ (alt 69°) to 1.8′ (alt 36°) of systematic offset after fixing B1. Drop the Python-side call and pass the true centre altitude; the kernels already handle the flag. (Inherited from the source branch.)

B3. GeoidHeightEGM96 signature change breaks nightly processing for cams_code > 0 stations. Four prerelease callers still pass a Config (CAL.py:119, Conversions.py:516,528,577,578, FrameInterface.py:1697). Verified: it raises AttributeError: 'Config' object has no attribute 'flush' from np.fromfile. writeCAL is called unguarded from Reprocess.py:636 whenever config.cams_code > 0, so processNight will raise for every CAMS-compatible station. This has to be fixed in this PR rather than left as a known gap. Simplest: accept None | str | Config (if hasattr(x, 'egm96_path'), build the path), which is what 7b5125e4 did before it was reverted.

B4. GPUENHt.py raises NameError, not ImportError, without numba. @cuda.jit at module level runs after the except ImportError handler, so cuda is undefined. Verified: import RMS.Astrometry.GPUENHt with numba absent → NameError: name 'cuda' is not defined; Janus's except ImportError will not catch it. Wrap the kernel in if CUDA_AVAILABLE:, or better, keep this pipeline-specific CUDA module out of RMS (no RMS caller, prints at import time).

Important

I1. Platepar.height_wgs84 is a persisted derived value that goes stale. Verified: Platepar() has no height_wgs84 (every ported transform raises AttributeError on a fresh platepar); pp.elev += 1000 leaves it unchanged, and Reprocess.py:120, ApplyRecalibrate.py:730,844 and SkyFit2's station dialog all overwrite elev; jsonStr() now writes it into every platepar file, and once present loadFromDict never recomputes it, so a stale value is baked in. Suggest not storing it: compute it inside the wrappers from lat/lon/elev, or make it a read-only @property stripped in jsonStr.

I2. mslToWGS84Height rebuilds the geoid spline on every call. interpolateEGM96Data (721×1440 RectSphereBivariateSpline) costs ~0.09 s here per distinct argument tuple; @memoizeSingle caches only the last call. Platepar.read() went from ~0 to 0.09 s on a desktop; expect roughly 0.5–1 s per load on a Pi 4. Cache the model per file path in loadEGM96Data/interpolateEGM96Data instead.

I3. cyXYHttoENU_wgs84 drops the poly3+radial3 / poly3+radial5 terms that cyXYToAltAz and the three *_iter kernels include, so XY→ENU and ENU→XY are not inverses for those distortion types. It also silently falls through to "no distortion" for an unrecognised dist_type.

I4. No tests. 4ab9bbbf removed the ENU tests and the branch's coordinate_transform_tests/ were not brought over. Given B1/B2, a test comparing xyToAltAzPP against xyToRaDecPP → alt/az on a shipped template platepar (plus a rolled variant) to < 0.5′ is the minimum.

I5. cyENHt0ToENHt1 / cyENHtToXY_iter bracketing. Bracket expansion moves U_lo but only re-evaluates f_hi, so f_lo goes stale. The second bisection in cyENHt0ToENHt1 assumes height increases monotonically along the ray and converges to t_hi = 1e6 m (garbage, not NaN) when Ht1 is never reached; a downward ray bisects in the wrong direction. Return NaN on bracket failure as cyXYHttoENU_wgs84 does.

Minor

• imageCenter and _centreAltAz are dead code here; _centreAltAz would also be wrong if used (it precesses epoch-of-date RA_d/dec_d).
• geodetic_to_ecef, ecef_to_geodetic_bowring, R_ecef_from_enu are unused.
• cyXYToAltAz lost the @cython.boundscheck/wraparound/cdivision decorators the source has; r_corr, r_scale (and dX, r_corr in cyXYHttoENU_wgs84) are undeclared and become Python objects in the hot loop. Its docstring says "az_centre: Reference right ascension".
• ENHt0ToENHt1 docstring promises "ndarray or float" for Ht1_data, but a scalar raises ValueError (verified); broadcast instead.
• xyHtToENUPP: np.full(len(X_data), …) fails if X_data is a scalar.
• enHtToXYPP imports cyENHtToXY_iter locally; move it to the module import block.
• xyToAltAzPP docstring documents a non-existent extinction_correction kwarg.
• GeoidHeightEGM96 CLI lost -c/--config.

What checks out

The seven kernels (apart from cyXYToAltAz's decorators) are byte-identical to the source branch; prerelease functions are untouched; XY→ENU→geodetic→XY round-trips to 0.015 px and the ellipsoid intersection hits the target height to 9 mm; memoizeSingle, getRmsRootDir and share/WW15MGH.DAC all exist on prerelease; Cython builds cleanly from a cold cache.

Happy to share the comparison scripts if useful.

[dvida]

B1 (blocking). rotationWrtHorizon returns atan2(dAlt, dAz) without cos(alt) on the azimuth step, so it is not the tangent-plane position angle the alt/az gnomonic in cyXYToAltAz needs. Measured max error vs xyToRaDecPP → alt/az: 28′–38′ on unrolled real platepars, up to 18.4° with pos_angle_ref + 45°. The best-fit rotation is exactly pos_angle_ref − q with q the parallactic angle at the FOV centre (atan2(sin HA, tan(lat)·cos(dec) − sin(dec)·cos(HA)), HA from JD2LST(pp.JD, pp.lon) and pp.RA_d), matching to < 0.001° on both test platepars; with that the max error drops to 0.1′–0.2′. Same issue at lines 1745, 1838, 1886, 1939.

[dvida]

B2 (blocking). This refracts the centre unconditionally, and cyXYToAltAz applies refractionTrueToApparent(h0) again when refraction=True (and applies nothing when False, so the centre is then refracted once where it should not be). Net effect: one extra refraction of the centre either way, 0.47′ at alt 69° to 1.8′ at alt 36°. Drop this line (and the same line in the other four wrappers: 1742, 1835, 1883, 1936) and pass the true altitude; the kernels already handle the flag.

[dvida]

B3 (blocking). Four prerelease callers still pass a Config here (CAL.py:119, Conversions.py:516,528,577,578, FrameInterface.py:1697). Verified: mslToWGS84Height(lat, lon, h, config) → AttributeError: 'Config' object has no attribute 'flush' from np.fromfile. writeCAL is called unguarded from Reprocess.py:636 for every station with cams_code > 0, so this breaks nightly processing for them. Please accept None | str | Config (e.g. if hasattr(egm96, 'egm96_path'): path = os.path.join(egm96.egm96_path, egm96.egm96_file_name)) rather than leaving it as a known gap.

[dvida]

I2. interpolateEGM96Data builds a 721×1440 RectSphereBivariateSpline on every call (~0.09 s on a desktop, likely 0.5–1 s on a Pi 4), and @memoizeSingle on the outer function only caches the last argument tuple, so every Platepar.read() of a different platepar pays the full cost. Cache the loaded array and the interpolated model per file path (memoize loadEGM96Data and interpolateEGM96Data) instead.

[dvida]

B4 (blocking). When numba is missing the except ImportError above sets CUDA_AVAILABLE = False, but this decorator still executes and raises NameError: name 'cuda' is not defined (verified). Janus's except ImportError guard will not catch a NameError. Wrap the kernel definition in if CUDA_AVAILABLE: or, preferably, keep this contrail-pipeline CUDA module outside RMS since nothing in RMS calls it.

[dvida]

I3. This poly3 branch omits the +radial3 / +radial5 extra terms (x_poly_fwd[12], [13]) that cyXYToAltAz, cyGeoToXY_wgs84_iter, cyENUToXY_iter and cyENHtToXY_iter all include, so XY→ENU and ENU→XY are not inverses for poly3+radial3 / poly3+radial5 platepars. The else at line 2187 also silently applies no distortion for an unknown dist_type instead of raising.

[dvida]

I5. When the bracket is expanded on the U_lo side only f_hi is re-evaluated, so f_lo is stale (same pattern at line 2927 in cyENHtToXY_iter). Further down, the second bisection (t_lo = 10, t_hi = 1e6) assumes height increases monotonically along the ray and, if Ht1 is never reached, converges to t_hi and returns that as a valid result rather than NaN; a downward ray bisects in the wrong direction. Please return NaN on bracket failure as cyXYHttoENU_wgs84 does.

[dvida]

Minor: the docstring says Ht1_data may be a float, but a scalar Ht1 raises this ValueError (verified: Got E0:3, N0:3, Ht0:3, Ht1:1). Broadcast scalars to the length of E0 instead.

[dvida]

Minor: _centreAltAz and imageCenter above have no callers on this branch. If _centreAltAz were used it would also be wrong: cyTrueRaDec2ApparentAltAz precesses J2000→date, but RA_d/dec_d are already epoch-of-date in RMS (see pointingCorrection), which is why the cyraDec2AltAz call in the wrappers is the right one.

[dvida]

Minor: the source branch has @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) on this function; they were dropped here, so it is not quite verbatim. r_corr and r_scale are also undeclared (Python objects in the loop), and the docstring describes az_centre/alt_centre as RA/Dec. The three cdef inline helpers below (geodetic_to_ecef, ecef_to_geodetic_bowring, R_ecef_from_enu) are unused.

[Cybis320]

A little background on rotationWrtHorizon (issue B.1): The angle has never been correct in master, and it did not really matter as the RaDec functions don't really care about the value being correct - they just require an anchor. However, the more efficient direct XY <> AltAz transforms copied here from my implementation do require an accurate value. My test-coordinate-transforms branch has the correct rotationWrtHorizon implementation and I don't believe it breaks current platepars.

Regarding issue B.2, it was fixed in test-coordinate-transforms if you want to port that over.

[liam-schultz]

If RMS is installed in a conda environment this is at os.path.join(sys.prefix, "share", "WW15MGH.DAC") instead. We should check both.