Fixes a north–south flip of the geopotential cube in the ECMWF/ERA-5 model-level reader.

[jlmaurer]

Description

Thanks to @s-sasaki-earthsea-wizard for catching this bug!

Fixes a north–south flip of the geopotential cube in the ECMWF/ERA-5 model-level reader.

In ECMWF._load_model_level, the branch that corrects ECMWF's descending-latitude ordering reversed the wrong axis of z:

if lats[0] > lats[1]:
    z    = z[::-1]        # z is (level, lat, lon) -> reverses LEVEL
    lnsp = lnsp[::-1]     # (lat, lon)             -> reverses lat
    t    = t[:, ::-1]     # (level, lat, lon)      -> reverses lat
    q    = q[:, ::-1]     # -> reverses lat

z needs both axes reversed: the latitude axis to stay aligned with t/q/lnsp, and the level axis because z[0] is subsequently passed to _calculategeoh as the surface geopotential. Reversing only axis 0 satisfied the second requirement while silently violating the first, leaving the height field mirrored in latitude against the meteorology.

The visible symptom is columns whose pressure profile is displaced vertically by however much terrain differs between a latitude and its mirror — e.g. a sea-level node reporting constant pressure from −500 m up to 1900 m, which is ~1900 m of fictitious sea-level-density air worth roughly +490 mm of spurious delay.

Affects ERA-5, ERA-5T and HRES on model levels, which is the default level type for those models. The pressure-level path is untouched.

When this was introduced

Introduced in f81225b ("Update for ERA-5 API changes", #751, merged 2025-09-05, released in 0.6.0) — so the regression window is 0.6.0 onward, not the whole history.

The statement z = z[::-1] is older than the bug and was correct when written. What changed is the dimensionality of z underneath it:

	before f81225b	after
_makeDataCubes	np.squeeze(block['z'].values)[0, ...] → 2D (lat, lon)	block['z'].values.squeeze() → 3D (level, lat, lon)
z[::-1]	reverses latitude ✓	reverses level ✗
call site	_calculategeoh(z, lnsp)	_calculategeoh(z[0, :], lnsp, t, q)

The ERA-5 API change meant CDS would only return z at the surface, so RAiDER had to compute the full geopotential cube itself. That moved the [0, ...] surface-extraction out of _makeDataCubes and into the call site as z[0, :], which made the level reversal newly load-bearing — it is what puts the surface at index 0 — at exactly the moment the latitude reversal silently stopped happening. One statement went from serving one purpose to serving a different one, and reads as deliberate either way.

Fixes #806

Why this was not caught earlier

The error is antisymmetric in latitude, so it very nearly cancels in a domain average while being severe at individual points. Over a 3-degree box in southern California the mean ZTD moves only 2360 → 2340 mm, but the spread across the domain collapses from physically impossible to correct:

diagnostic	before	after
corr(surface height, −surface pressure)	−0.5415	+0.9992
sea-level ZTD range across the domain	1621 – 3233 mm	2239 – 2395 mm
domain-mean ZTD	2360 mm	2340 mm

Physics requires the correlation to be ~+1 (high terrain carries low surface pressure).

Validated against UNR GNSS ZTD at four stations (FGNW, MHMS, WLHG, WORG) on four dates in 2020 and 2025, spanning both dry and high-water-vapour conditions. Per-station residuals go from ±400–500 mm to −13 to −66 mm.

Tests

test/test_ecmwf_levels.py builds a synthetic ERA-5 model-level file in ECMWF's native descending-latitude ordering, with terrain rising from 0 m in the south to 3000 m in the north so a reflection is unambiguous, and loads it through the real _load_model_level. z is generated with calcgeoh exactly as RAiDER's own _fetch writes it.

Two assertions are load-bearing and both fail on the pre-fix code (the correlation comes out at −0.9986):

• surface height per latitude row matches the terrain it was built from
• corr(surface height, −surface pressure) > 0.99

The domain mean is deliberately not asserted on, for the reason above; every assertion is per-column.

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.

[s-sasaki-earthsea-wizard]

@jlmaurer Confirmed — this fixes the latitude mirror. I re-ran the reproduction harness from the evidence repo against this branch at its current head 919d6b6, with the same raw ERA-5 ML file and the same frozen dependency set as the dev@1fb9e14 and #805e7c8187 runs reported in #806 — only the RAiDER checkout changed (raider==0.6.1.dev136+g919d6b699).

Result: the mirror is gone.

diagnostic (575 columns, 4.5°×4.5° box over Kanto)	dev@1fb9e14 (= #805e7c8187)	this PR @919d6b6
z_surface passed to _calculategeoh (bitwise ID)	lowest full level, latitude-mirrored	lowest full level, own latitude
height offset vs true-surface reference	−1252.41 / +1271.91 m	+9.46 / +10.19 m
regression dh = a + b·h_own + c·h_mirror	b = −1.000243, c = +0.999979	b = −0.000243, c = −0.000021

The fitted coefficients changed by exactly +1.000000 for h_own and −1.000000 for h_mirror, removing the mirror component (−h_own + h_mirror) while the residual component (intercept 9.867 m, rmse 0.149 m) is unchanged to the last digit, and the measured post-fix offset lands inside the minimal-fix simulation from the #806 report (predicted +9.44 to +10.20 m). This matches the simulated axis-only fix quantitatively.

I also ran the new test/test_ecmwf_levels.py in both environments: 4/4 pass on this branch, and on pre-fix dev@1fb9e14 the two load-bearing assertions fail (correlation comes out at −0.9986, matching the PR description). The descending-latitude fixture closes the test-suite gap I noted on #805 — this regression can't silently return now.

Updated results are in the evidence repo, pinned at 67087c0: see reports/ml_offset_regression_pr811.md.

The remaining ~10 m offset (tracked separately)

The remaining ~+9.5–10.2 m is a separate, much smaller correctness issue: after the fix, z[0] is the lowest full model level where _calculategeoh expects the true surface field, so the reconstructed columns retain a near-uniform vertical offset of that size (measured above). This is the "additional item" from my #806 proposal and I don't think it should block this focused fix — it is about two orders of magnitude below the bug fixed here and roughly uniform in space. Since merging this PR closes #806, I've opened #812 to track the separate z_surface + loader-fallback change, including how it interacts with the reworked fetch in #805. Happy to follow up with that PR once this merges.

[sssangha]

Great catch on this and thanks @s-sasaki-earthsea-wizard for the detailed verification!

Quick review notes where I distilled everything for easy reference later:

• The Root Cause: Spot-on fix. When z became a 3D cube in Fix for ERA-5 API change #751, z[::-1] started quietly flipping the vertical level axis instead of latitude, leaving z mirrored North-South against t, q, and lnsp.
• Validation: The diagnostics here look good. Seeing UNR GNSS ZTD residuals drop from ±400–500 mm down to -13 to -66 mm, and the surface terrain/pressure correlation jump back to +0.9992 shows the distortion is addressed.
• Testing: Good we captured the flip in test/test_ecmwf_levels.py and made it clear in the notes.

LGTM! I'll approve

[sssangha]

LGTM, see my comment

[s-sasaki-earthsea-wizard]

Thanks @sssangha !
That summary is a nice reference to have in one place.