Summary
GdalRaster::createTransform never sets a coordinate epoch, so when the dataframe CRS is dynamic (ITRF2020 EPSG:9989 for r007 products via #534; the EPSG:7912 ITRF2014 fallback likewise) and the sampled raster's CRS is static (e.g. 3DEP 1 m: NAD83(2011) UTM + NAVD88), PROJ evaluates the time-dependent Helmert ("ITRF2020 to NAD83(2011) (1)", EPSG:10334, t_epoch=2010.0) at its parameter reference epoch 2010.0 instead of the observation epoch. The sampled location and the vertical-shift correction are therefore both computed for where the North America plate was in 2010.
Measured effect (western US, pyproj 3.7.2 / PROJ 9.7.1)
Difference between the 4D transform EPSG:9989 → EPSG:6318+5703 evaluated at the observation epoch vs at 2010.0, on a 362-point grid across the western US (exactly linear in time; azimuth is the direction the sampled location is displaced from the epoch-correct one):
| epoch |
horizontal error (m) |
vertical error of the POI orthometric height (mm) |
| 2018.75 |
0.169 |
+9.2 |
| 2022.75 |
0.246 |
+13.4 |
| 2026.60 |
0.320 |
+17.5 |
Rates: 19.3 mm/yr horizontal toward azimuth ~238° (rigid NA-plate motion; varies ~1–2.5 cm/yr with location across CONUS since it is a plate rotation) and +1.05 mm/yr vertical. Because samplePOI computes vertical_shift = z_source − z_transformed (GdalRaster.cpp:331) and adds it to the pixel value, the vertical part propagates directly into the returned value: 3dep.value is biased low by ~1.5 cm for 2024 data, growing ~1 mm/yr — a spurious trend for anyone using SlideRule for multi-year elevation change against 3DEP. The horizontal part samples the DEM ~0.2–0.3 m away from the correct pixel: on terrain slope s this adds a per-point error ≈ displacement × s (for our mountain AOI: ~3 cm NMAD-equivalent scatter at median slopes, 0.15–0.2 m per-point on slopes > 0.5 by 2026).
Verification context: we verified the rest of the chain is exact — client-side resampling of the same 3DEP COGs with an epoch-aware 4D PROJ transform (GEOID18 + the same Helmert) reproduces 3dep.value to 0.000 m median/NMAD once the epoch is forced to 2010.0. May be relevant background for other standing-offset reports against 3DEP (e.g. #263).
Code pointers (current main, apps/node/ layout)
apps/node/packages/geo/package/GdalRaster.cpp:1286–1358 (createTransform): OGRCoordinateTransformationOptions gets an AOI / optional user pipeline, but no coordinate epoch is ever set before OGRCreateCoordinateTransformation(&sourceCRS, &targetCRS, options).apps/node/packages/geo/package/GdalRaster.cpp:304–368 (samplePOI): transforms the POI and derives vertical_shift (line 331), added to the sampled value (e.g. lines 902, 974).apps/node/packages/geo/package/RasterObject.h:75–79: point_info_t already carries a per-point int64_t gps timestamp — the epoch information is plumbed to the sampling layer.
Suggested fix
Call sourceCRS.SetCoordinateEpoch(epochYears) (GDAL ≥ 3.4, RFC 81 — the coordinate epoch is a property of OGRSpatialReference, not of OGRCoordinateTransformationOptions) in createTransform before creating the transformation, with epochYears = decimal year of the request's mid-time (e.g. from the first/median point_info_t.gps of the batch, or granule metadata). GDAL then instantiates the PROJ operation at that epoch.
Granularity: within a granule-region request the point times span seconds, so one epoch per transform is exact for practical purposes; even 0.1-yr rounding (if useful for transform caching) keeps the residual under ~2 mm horizontal. Transforms reused across requests (if any) would need the epoch in their cache key. Static-source dataframes are unaffected: setting an epoch on a static CRS is a no-op for static-to-static operations, though guarding with sourceCRS.IsDynamic() keeps intent explicit.
Workaround for users (until fixed)
The vertical component can be removed post-hoc: value_corrected = value + 1.05 mm/yr × (t_obs − 2010.0) for CONUS (location-dependent at the ±10% level). The horizontal (slope-coupled) component cannot be repaired after sampling.
Full analysis scripts available on request (pyproj 4D grid evaluation + COG-resampling verification).
Summary
GdalRaster::createTransformnever sets a coordinate epoch, so when the dataframe CRS is dynamic (ITRF2020EPSG:9989for r007 products via #534; theEPSG:7912ITRF2014 fallback likewise) and the sampled raster's CRS is static (e.g. 3DEP 1 m: NAD83(2011) UTM + NAVD88), PROJ evaluates the time-dependent Helmert ("ITRF2020 to NAD83(2011) (1)", EPSG:10334,t_epoch=2010.0) at its parameter reference epoch 2010.0 instead of the observation epoch. The sampled location and the vertical-shift correction are therefore both computed for where the North America plate was in 2010.Measured effect (western US, pyproj 3.7.2 / PROJ 9.7.1)
Difference between the 4D transform
EPSG:9989 → EPSG:6318+5703evaluated at the observation epoch vs at 2010.0, on a 362-point grid across the western US (exactly linear in time; azimuth is the direction the sampled location is displaced from the epoch-correct one):Rates: 19.3 mm/yr horizontal toward azimuth ~238° (rigid NA-plate motion; varies ~1–2.5 cm/yr with location across CONUS since it is a plate rotation) and +1.05 mm/yr vertical. Because
samplePOIcomputesvertical_shift = z_source − z_transformed(GdalRaster.cpp:331) and adds it to the pixel value, the vertical part propagates directly into the returned value:3dep.valueis biased low by ~1.5 cm for 2024 data, growing ~1 mm/yr — a spurious trend for anyone using SlideRule for multi-year elevation change against 3DEP. The horizontal part samples the DEM ~0.2–0.3 m away from the correct pixel: on terrain slopesthis adds a per-point error≈ displacement × s(for our mountain AOI: ~3 cm NMAD-equivalent scatter at median slopes, 0.15–0.2 m per-point on slopes > 0.5 by 2026).Verification context: we verified the rest of the chain is exact — client-side resampling of the same 3DEP COGs with an epoch-aware 4D PROJ transform (GEOID18 + the same Helmert) reproduces
3dep.valueto 0.000 m median/NMAD once the epoch is forced to 2010.0. May be relevant background for other standing-offset reports against 3DEP (e.g. #263).Code pointers (current main,
apps/node/layout)apps/node/packages/geo/package/GdalRaster.cpp:1286–1358(createTransform):OGRCoordinateTransformationOptionsgets an AOI / optional user pipeline, but no coordinate epoch is ever set beforeOGRCreateCoordinateTransformation(&sourceCRS, &targetCRS, options).apps/node/packages/geo/package/GdalRaster.cpp:304–368(samplePOI): transforms the POI and derivesvertical_shift(line 331), added to the sampled value (e.g. lines 902, 974).apps/node/packages/geo/package/RasterObject.h:75–79:point_info_talready carries a per-pointint64_t gpstimestamp — the epoch information is plumbed to the sampling layer.Suggested fix
Call
sourceCRS.SetCoordinateEpoch(epochYears)(GDAL ≥ 3.4, RFC 81 — the coordinate epoch is a property ofOGRSpatialReference, not ofOGRCoordinateTransformationOptions) increateTransformbefore creating the transformation, withepochYears= decimal year of the request's mid-time (e.g. from the first/medianpoint_info_t.gpsof the batch, or granule metadata). GDAL then instantiates the PROJ operation at that epoch.Granularity: within a granule-region request the point times span seconds, so one epoch per transform is exact for practical purposes; even 0.1-yr rounding (if useful for transform caching) keeps the residual under ~2 mm horizontal. Transforms reused across requests (if any) would need the epoch in their cache key. Static-source dataframes are unaffected: setting an epoch on a static CRS is a no-op for static-to-static operations, though guarding with
sourceCRS.IsDynamic()keeps intent explicit.Workaround for users (until fixed)
The vertical component can be removed post-hoc:
value_corrected = value + 1.05 mm/yr × (t_obs − 2010.0)for CONUS (location-dependent at the ±10% level). The horizontal (slope-coupled) component cannot be repaired after sampling.Full analysis scripts available on request (pyproj 4D grid evaluation + COG-resampling verification).