Skip to content

Fix ellipsoidal/geoid height datum handling for station and DEM queries - #822

Open
royagrace wants to merge 9 commits into
dbekaert:devfrom
royagrace:stations
Open

Fix ellipsoidal/geoid height datum handling for station and DEM queries#822
royagrace wants to merge 9 commits into
dbekaert:devfrom
royagrace:stations

Conversation

@royagrace

@royagrace royagrace commented Sep 2, 2026

Copy link
Copy Markdown
Contributor
  1. dem.py / llreader.py: DEM heights and GNSS station heights are now correctly normalized to geoid (~MSL) to match ERA5's convention, fixing the original ZTD sampling bias.
  2. llreader.py / delay.py / losreader.py: readZ(ellipsoidal_heights=...) now lets callers get either geoid (for weather-model sampling) or true ellipsoidal (for LOS/ECEF geometry) heights, and delay.py requests the right one for each consumer instead of reusing one value for both.
  3. A few stale "ellipsoidal" comments/docstrings elsewhere (utilFcns.py, weatherModel.py, delay.py) were corrected to say geoid, matching what the code actually does.
  4. Test coverage added throughout for the new conversion paths.

Description

  • dem.py: download_dem now requests geoid-referenced GLO-30 DEM heights (dst_ellipsoidal_height=False) instead of ellipsoidal, matching ERA5's convention.
  • llreader.py:
    • Added _ellipsoidal_to_geometric/_geometric_to_ellipsoidal, a bidirectional height-datum converter (via EPSG:9707, WGS84+EGM96), with a PROJ-CDN fallback if the local EGM96 grid is missing.
    • StationFile now accepts a crs param (default 4326=geoid; pass 4979 for ellipsoidal GNSS heights, e.g. UNR MAGNET/IGS20) and update_crs() to change it after construction.
    • readZ() on StationFile, RasterRDR, and GeocodedFile takes a new ellipsoidal_heights flag (default False=geoid) so callers can request either convention; each converts from its known source datum on demand.
  • delay.py: tropo_delay now reads heights twice with the correct flag for each use — geoid for transformPoints/weather-model sampling, ellipsoidal for los.setPoints/LOS-ECEF geometry — instead of reusing one value for both, which was silently wrong for whichever consumer didn't get its expected convention.
  • losreader.py: LOS.setPoints() renamed heightsellipsoidal_heights and documented that it requires true WGS84-ellipsoidal input, since it feeds state_to_los/lla2ecef, which is geodetic-to-ECEF math and only correct with ellipsoidal height.
  • Stale comments fixed (no behavior change) in utilFcns.py (geo_to_ht), models/weatherModel.py (_get_heights), and delay.py (cube z attrs, transformPoints docstring) — these compute/label geoid-referenced heights but previously said "ellipsoidal."
  • Test coverage added across all of the above (test_llreader.py).

Needs revisiting

  • RasterRDR/GeocodedFile height-source datum is assumed, not verified. When heights come from a user-supplied hgt_file or an existing is_dem file (rather than a freshly-downloaded DEM), there's no datum actually tracked — we assume geoid to match everything else and documented that assumption inline, but it hasn't been confirmed for real-world files of that kind.
  • self._crs ellipsoidal/geoid detection is a heuristic, not a real check: StationFile infers "ellipsoidal" from len(crs.axis_info) >= 3 (does the CRS have a vertical axis), not from the CRS's actual datum. This is correct for the two CRSs the API documents (4326, 4979), but would misclassify a 3D geoid CRS like EPSG:9707 as ellipsoidal if someone passed that in.
  • RasterRDR.readZ(ellipsoidal_heights=True) raises NotImplementedError if the object was built without a lon_file (the incomplete 2-band-raster case), rather than supporting it — flagged, not fixed.
  • Cube/raytracing path (BoundingBox/Geocube, _build_cube_ray) not audited. This PR only fixes the point-AOI path (StationFile/RasterRDR/GeocodedFile) that goes through delay.py's los.setPoints(). Whether the raytracing cube-generation path has the same ellipsoidal/geoid mismatch wasn't checked.

Motivation and Context

RAiDER was sampling the weather-model delay cube with GNSS station heights that were ellipsoidal (WGS84), while the cube's own z-axis is geoid-referenced (~MSL), like ERA5. This PR normalizes height handling so every height is explicitly geoid or ellipsoidal, and each consumer gets the one it actually needs.
Fixes some issues from PR #805 and all issues from PR #813

How Has This Been Tested?

Ran autopep8 command
Created and ran tests locally in test_llreader.py

Screenshots (if appropriate):

Type of change

  • Bug fix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)
  • Breaking change (fix or feature that would cause existing functionality to not work as expected)

Checklist:

  • I have added an explanation of what your changes do and why you'd like us to include them.
  • I have written new tests for your core changes, as applicable.
  • I have successfully ran tests with your changes locally.
  • My change requires a change to the documentation.
  • I have updated the documentation accordingly.

@royagrace
royagrace marked this pull request as ready for review September 2, 2026 20:18
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant