VDatum mosaic: offsets of a few cm where packages from different generations meet (visible in vdatum:msl) #192

I found this while testing the fix for #190 (PR #191). This bug was already there before that work; it may have come from the VDatum mosaicking added back in #161.

What's happening

Transforms through vdatum:msl near the edge of a VDatum package have patches of cells offset by up to ~5 cm.

Repro: region [-76.5, -76.4, 38.9, 39.0] (Annapolis area), 3 arc-second grid.

• mtl -> msl has a median of -0.006 m, but 517 of 14,280 cells--in several patches along the North and West edges over land--make sudden jumps to -0.026 to -0.056 m.
• mllw -> msl shows a similar, slightly larger patch of 616 cells.
• mtl -> mllw and mtl -> navd88 do not.

To build the grids above and inspect them yourself:

transformez build -R -76.5/-76.4/38.9/39.0 -E 3s -I vdatum:mtl -O vdatum:msl -o mtl_to_msl.tif
transformez build -R -76.5/-76.4/38.9/39.0 -E 3s -I vdatum:mllw -O vdatum:msl -o mllw_to_msl.tif
transformez build -R -76.5/-76.4/38.9/39.0 -E 3s -I vdatum:mtl -O vdatum:mllw -o mtl_to_mllw.tif
transformez build -R -76.5/-76.4/38.9/39.0 -E 3s -I vdatum:mtl -O epsg:5703 -o mtl_to_navd88.tif

Cause

Three packages cover this region, in priority order:

1. MDnwchb11_8301: 2024 release, xGEOID20B/IGS14 roadmap.

2. MDVAechb11_8301: TSS is empty here.

3. MDVAchb12_8301: 2019 release, NAVD88 roadmap.

4. MSL is built from a different set of cells than the other tidal surfaces. _fetch_vdatum_coverage_paths builds MSL from the TSS grids alone. Every other surface needs both the tidal grid and the TSS grid from the same package. TSS grids extend further than the tidal grids, mostly onto land. For MSL, MDVAchb12 therefore fills 102 cells that MDnwchb11 lacks. For mtl it contributes 0 cells, so those cells get the global:mss blend instead. coastal_aware_composite then spreads the offset into ~400 neighbouring cells.

5. The two packages don't agree with each other. Where both packages have data, MDVAchb12 sits above MDnwchb11 after normalization, for every surface:

   surface	mean	max
   msl	+2.5 cm	4.9 cm
   mtl	+2.6 cm	5.3 cm
   mllw	+3.5 cm	6.8 cm
   mhw	+2.0 cm	4.9 cm

   It breaks down into two parts:

   • ~2.5 cm shared by all surfaces. VDatum's own TSS grids put NAVD88 0.386 m above xGEOID20B here. The normalization chain (xGEOID20B → IGS14@2010 → NAD83(2011) via HTDP → minus GEOID18) gives 0.361 m.
   • A part that differs by surface, from the packages' own tidal models: mtl +0.2 cm, mllw +1.2 cm, mhw −0.4 cm.

   The Annapolis tide gauge (8575512, inside the region) suggests this is a disagreement in the source data, not a normalization bug:

   NAVD88 − MSL	value
   CO-OPS published	+0.016 m
   MDVAchb12	+0.019 m
   MDnwchb11 via the normalization chain	+0.065 m

   For MLLW − MSL (gauge −0.219 m), the newer package is closer: MDnwchb11 gives −0.221 m, MDVAchb12 gives −0.209 m. Each generation is consistent within its own system, and GEOID18 doesn't reconcile them at this location.

Ruled out:

• Package selection or ordering: all surfaces get the same packages in the same order.
• Inland decay: extrapolate_inland=True gives identical results.
• xGEOID20B epoch: vdatum:xgeoid20b uses EPSG:7912 at epoch 2010, while NGS publishes xGEOID20B at 2020.0. Using 2020 moves the result 0.8 cm further from the gauge. Also, HTDP 3.6.0 reports these points as "outside of the modeled region" whenever the epoch changes. The epoch is probably worth correcting on its own, but it isn't the cause here.

MSL is where this showed up in my tests. The same seam will appear for any tidal surface wherever a lower-priority package from another generation fills cells.

Possible fixes

A. Build MSL only from cells where the same package also has tidal grids.

• Pros: small change; makes MSL consistent with the other surfaces; fixes this repro (dropping MDVAchb12 from MSL removes the outliers).
• Cons: MSL loses VDatum coverage on those edge cells, which fall back to the global blend. It doesn't help other surfaces where packages from different generations meet.

B. Shift a lower-priority package by its mean offset from the higher-priority package before letting it fill cells. Measure the offset where both have data, ideally near the cells being filled (2.9 cm here).

• Pros: general, covering all surfaces and generation mixes; keeps VDatum coverage; removes most of the seam.
• Cons: a single constant can't remove the part that varies across the area (the gap ranges from ~4 cm NW to ~1.4 cm S). It moves the lower-priority package away from its own gauge ties. It needs a minimum overlap and a fallback when there is none.

C. Don't let a package from a different generation fill cells next to a higher-priority one. Those cells would use the global/decay fill instead.

• Pros: simple; no mixed-generation seams at all.
• Cons: discards valid VDatum data. The global fill may be less accurate than the package being dropped, especially in estuaries.

D. Blend the overlap region instead of switching sharply from one package to the other.

• Pros: turns the step into a smooth change and keeps all the data.
• Cons: the blend zone is correct under neither package. It adds complexity to _build_vdatum_hydro_mosaic, and it doesn't help where there's no overlap.

Suggestion: A now, since it's low-risk and fixes the visible MSL symptom, then choose between B, C, or D (or some other solution) after that.

I'm not sure the best approach here @matth-love, but I wanted to bring it up for your review and to find a possible fix.

Also noted (but probably unrelated):

The 4 commands above each produce transformation grids with 120 rows (3" spacing) and 119 columns (3.025" spacing), even though 3s grids over 0.10 degrees would be expected to contain 120 cells of exactly 3" in both dimensions.