Systematic offset in elevations returned between GEDI and 3DEP

[jpswinski]

@dshean: "[There is] a systematic offset in the elevations returned for GEDI L4A and 3DEP, which indicates there is a residual issue with the vertical datum handling, or our assumption about GEDI product vertical datum is incorrect. An offset of ~15 m is consistent with the NAVD88/geoid height and the height above ellipsoid for Grand Mesa."

"This ties in with all of the discussion around the PROJ transformation and the geoid_grid offset we’ve been discussing for 3DEP support. If we are doing the transformation properly, the difference values should center around 0."

"A systematic offset of 15 m is a vertical datum problem. That’s the magnitude of the difference between the geoid and the ellipsoid height for Grand Mesa."

[jpswinski]

[jpswinski]

We are likely improperly handling the CRS of one of the two datasets, or we are improperly specifying the transformation required to convert one to the other. Needs further investigation.

[jpswinski]

@dshean: "I see the same offset for ICESat-2 returns limited to ATL08 ground photons. Looks like we’re not accounting for the vertical datum properly. Perhaps an issue accessing the CDN geoid offset grids?"

[scottyhq]

Approach 1:

1. Warp the raster (via VRT), 2. then sample with points in EPSG:7912

# Virtual Warp raster to 7912 
AWS_NO_SIGN_REQUEST=YES gdalwarp -s_srs epsg:6342+5703 -t_srs EPSG:7912 /vsis3/prd-tnm/StagedProducts/Elevation/1m/Projects/CO_MesaCo_QL2_UTM12_2016/TIFF/USGS_one_meter_x23y434_CO_MesaCo_QL2_UTM12_2016.tif warped.vrt

# Sample (will actually apply transform and grid shifts behind the scenes)
echo -108.1 39.1 | CPL_DEBUG=ON PROJ_NETWORK=ON PROJ_DEBUG=2 AWS_NO_SIGN_REQUEST=YES GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR gdallocationinfo -geoloc warped.vrt

Report:
  Location: (2414P,6744L)
Value: 2638.0322265625

Approach 2:

1. Transform the IS2 Points to Raster CRS, 2. sample, 3. manually apply grid shift. (Lon_i, Lat_i, Z_i, t_i) -> ( Y_i, X_i, Z_i, t_i) where elevation Zi and time of observation t_i are optional

# 1. Transform 7912 to raster CRS 
# PROJ (cs2cs) or GDAL/OGR for the 3D transform 
# NOTE: Not sure how much in practice just using 0.0 for all elevations affects this 3D transform (likely depends on transform and location of input point)...
echo -108.1 39.1 2630.0 | PROJ_DEBUG=2 PROJ_NETWORK=ON cs2cs -f "%.3f" -r epsg:7912 epsg:6342+5703
231917.919  4332449.531 2646.526

# So Geoid shift is (Z0-Z1) = -16.526

# 2. Sample with raster with transformed coords (in this case UTM)
echo 231917.92 4332449.53   | CPL_DEBUG=ON PROJ_NETWORK=ON PROJ_DEBUG=2 AWS_NO_SIGN_REQUEST=YES GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR gdallocationinfo -geoloc /vsis3/prd-tnm/StagedProducts/Elevation/1m/Projects/CO_MesaCo_QL2_UTM12_2016/TIFF/USGS_one_meter_x23y434_CO_MesaCo_QL2_UTM12_2016.tif 
Location: (1923P,7556L)
Value: 2654.55810546875

# 3. should equal result of approach 1, but possible rounding and performance differences 
RETURN_VAL = 2654.55810546875 + -16.526 =  2638.032  

[dshean]

Thanks @scottyhq! At some point, it would be useful to benchmark these two approaches based on point sample sizes and count of raster tiles. My intuition is that approach #2 will be much more efficient for most typical use cases (AOI size, number of tiles, number of points).

[dshean]

NOTE: Not sure how much in practice just using 0.0 for all elevations affects this 3D transform (likely depends on transform and location of input point)...

If the point is close to ellipsoid surface, this is negligible. But for points several km above ellipsoid surface, the output horizontal position will be slightly different. Highest elevations in CONUS (for this 3DEP specific transform case) are ~4.4 km, but other cases could be >7-8 km (e.g., Himalayas). We should just use the 3D coordinates by default.

Example for Mt. Whitney, highest point in CONUS

echo -118.29 36.58 0 | PROJ_DEBUG=2 PROJ_NETWORK=ON cs2cs -f "%.3f" -r epsg:7912 epsg:6342+5703 -692023.806 4131532.366 25.875

echo -118.29 36.58 4421 | PROJ_DEBUG=2 PROJ_NETWORK=ON cs2cs -f "%.3f" -r epsg:7912 epsg:6342+5703 -692023.808 4131532.367 4446.875

In this case, coords are different by a few mm, but some other cases I tested at different CONUS locations were ~1 cm.

[jpswinski]

The GEDI to 3DEP elevations compare very well; as do the comparison between ICESat-2 and ArcticDEM mosaic. But there is still an offset between the ICESat-2 elevations and the ArcticDEM strips. Investigating as to the root cause.

[dshean]

Great! Let's make sure we are doing the assessment using ICESat-2 points over exposed, stable ground surfaces, not changing ice or vegetation.
And consider the median vertical offset (IS2 - sampled ArcticDEM/REMA strip values) for a large number (>10) of strips at the same location, not just a single strip. Individual strips can have absolute vertical error of +/-5 m, but this should be random, with no bias (in principle).