3DEP Seamless DEMs - 3D Coordinate Reference System Transformation Resources

Horizontal Coordinates: NAD83 [EPSG: 4269]
Vertical Coordinates: NAVD88 [EPSG: 5703]

https://portal.opentopography.org/raster?opentopoID=OTNED.012021.4269.1

More important references:

%%bash

CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 AWS_NO_SIGN_REQUEST=YES \
 gdalinfo -nofl /vsis3/prd-tnm/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt

Driver: VRT/Virtual Raster
Size is 3888006, 939612
Coordinate System is:
GEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101004,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4269]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000555556093587,72.000555556294898)
Pixel Size = (0.000092592592660,-0.000092592592660)
Corner Coordinates:
Upper Left  (-180.0005556,  72.0005556) (180d 0' 2.00"W, 72d 0' 2.00"N)
Lower Left  (-180.0005556, -15.0005556) (180d 0' 2.00"W, 15d 0' 2.00"S)
Upper Right ( 180.0000003,  72.0005556) (180d 0' 0.00"E, 72d 0' 2.00"N)
Lower Right ( 180.0000003, -15.0005556) (180d 0' 0.00"E, 15d 0' 2.00"S)
Center      (  -0.0002776,  28.5000000) (  0d 0' 1.00"W, 28d30' 0.00"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  NoData Value=-999999
  Overviews: 1944003x469806, 972001x234903, 486000x117451, 243000x58725, 121500x29362

%%bash
# NOTE: -r bilinear or cubic will report fractional pixel position
INPUT=/vsis3/prd-tnm/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 AWS_NO_SIGN_REQUEST=YES \
 gdallocationinfo -geoloc $INPUT -106.5000000 38.5000000

Report:
  Location: (793805P,361805L)
  Band 1:
    <LocationInfo><File>/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n39w107/USGS_13_n39w107.tif</File></LocationInfo>
    Value: 2759.5849609375

# NOTE: that 'EPSG:4326' and EPSG:4269 are both lon,lat but not equivalent!
# This is fairly uncertain, it's going from EPSG ensemble to NAD83(1986)
# Implicitly the 'top' choice of change is to do nothing and assume +/-4m acccuracy
#!projinfo -s EPSG:4326 -t EPSG:4269 -o proj --spatial-test intersects

3DEP Seamless metadata¶

Metadata files for this particular 10m tile live on AWS S3:

https://prd-tnm.s3.amazonaws.com/index.html?prefix=StagedProducts/Elevation/13/TIFF/current/n39w107/

# Look at metadata
import geopandas as gpd
gf = gpd.read_file('https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/TIFF/current/n39w107/USGS_13_n39w107.gpkg')

# Just treat polygon boundaries as EPSG:4326
gf = gf.to_crs('EPSG:4326')
gf

# Interesting! Going into this tile we have 13 different data sources. Some going back to 1958
# Which dataset is out query point accessing?
point=gpd.GeoSeries(gpd.points_from_xy(x=[-106.50000], y=[38.50000], crs='EPSG:4326'))
m = gf.explore(column='workpackage')
point.explore(m=m, color='magenta')

PROJ Transforms¶

# To WGS84 1150
# top transform is to just apply a vertical shift grid us_noaa_g2018u0.tif
!projinfo -s EPSG:6318+5703 -t EPSG:7661 -o PROJ --hide-ballpark --spatial-test intersects | head -n 20

Candidate operations found: 22
-------------------------------------
Operation No. 1:

unknown id, Inverse of NAD83(2011) to NAVD88 height (3) + NAD83(2011) to WGS 84 (1) + WGS 84 to WGS 84 (G1150), 4.015 m, United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 2:

unknown id, Inverse of NAD83(NSRS2007) to NAD83(2011) (1) + Inverse of NAD83(NSRS2007) to NAVD88 height (1) + NAD83(NSRS2007) to WGS 84 (1) + WGS 84 to WGS 84 (G1150), 4.1 m, United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming.

PROJ string:

# NOTE: override NAD83 (EPSG:4269) with NAD83(2011) (EPSG:6318)
!projinfo -s EPSG:6318+5703 -t EPSG:7912 -o PROJ --hide-ballpark --spatial-test intersects | head -n 30

Candidate operations found: 42
-------------------------------------
Operation No. 1:

unknown id, Inverse of NAD83(2011) to NAVD88 height (3) + Conversion from NAD83(2011) (geog3D) to NAD83(2011) (geocentric) + Inverse of ITRF2014 to NAD83(2011) (1) + Conversion from ITRF2014 (geocentric) to ITRF2014 (geog3D), 0.015 m, United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1
  +step +proj=cart +ellps=GRS80
  +step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138
        +ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006
        +dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05
        +ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame
  +step +inv +proj=cart +ellps=GRS80
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 2:

unknown id, Inverse of NAD83(NSRS2007) to NAD83(2011) (1) + Inverse of NAD83(NSRS2007) to NAVD88 height (1) + Null geographic offset from NAD83(NSRS2007) (geog3D) to NAD83(NSRS2007) (geog2D) + NAD83(NSRS2007) to NAD83(2011) (1) + Conversion from NAD83(2011) (geog2D) to NAD83(2011) (geocentric) + Inverse of ITRF2014 to NAD83(2011) (1) + Conversion from ITRF2014 (geocentric) to ITRF2014 (geog3D), 0.15 m, United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=gridshift +no_z_transform

Annotated pipeline¶

+proj=pipeline
  # Inverse of NAD83(2011) to NAVD88 height
  # (NAVD88 Geoid Height -> NAD83(2011) Ellipsoid Height)
  +step +proj=axisswap +order=2,1 
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 
  # Conversion from NAD83(2011) (geog3D) to NAD83(2011) (geocentric)
  +step +proj=cart +ellps=GRS80 
  # Inverse of ITRF2014 to NAD83(2011)
  # (NAD83(2011) 3D -> ITRF2014 3D)
  +step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138 
        +ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006
        +dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05
        +ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame 
  # Conversion from ITRF2014 (geocentric) to ITRF2014 (geog3D)
  +step +inv +proj=cart +ellps=GRS80
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

# Top ranked transform as single line to pass to GDAL
!projinfo -s EPSG:6318+5703 -t EPSG:7912 -o PROJ --hide-ballpark --spatial-test intersects --single-line -q

+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138 +ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006 +dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05 +ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

%%bash
# VRT of a VRT inception

SSRS='EPSG:9055+5773'
TSRS='EPSG:7912'
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
PROJ_PIPELINE='+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138 +ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006 +dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05 +ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1'
INPUT=/vsis3/prd-tnm/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt
OUTPUT=USGS_Seamless_DEM_13_7912.vrt

echo $PROJ_PIPELINE

CPL_DEBUG=OFF PROJ_DEBUG=2 \
 PROJ_NETWORK=ON \
 GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 AWS_NO_SIGN_REQUEST=YES \
 gdalwarp -overwrite -wm 500 -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 -te -180.0005556 -15.0005556 180.0000003 72.0005556 -tr 0.000092592592660 0.000092592592660 -r bilinear -ot Float32 -s_srs $SSRS -t_srs $TSRS -ct "${PROJ_PIPELINE}" ${INPUT} ${OUTPUT}

+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138 +ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006 +dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05 +ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
Creating output file that is 3888006P x 939612L.
Using internal nodata values (e.g. -999999) for image /vsis3/prd-tnm/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt.
Copying nodata values from source /vsis3/prd-tnm/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt to destination USGS_Seamless_DEM_13_7912.vrt.
Processing /vsis3/prd-tnm/StagedProducts/Elevation/13/TIFF/USGS_Seamless_DEM_13.vrt [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

%%bash

INPUT=USGS_Seamless_DEM_13_7912.vrt

# Value: 2744.60235816757

CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 PROJ_NETWORK=ON \
 PROJ_DEBUG=2 \
 AWS_NO_SIGN_REQUEST=YES \
 gdallocationinfo -r bilinear -geoloc $INPUT -106.5000000 38.5000000

Report:
  Location: (793805.999902109P,361806.000216605L)
  Band 1:
    Value: 2744.60235816757

%%bash
# Hmmmm, stated accuracy of transform is 0.015 m , but seeing a difference of 0.12m
PROJ_PIPELINE_7912='+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138 +ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006 +dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05 +ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1'
echo 38.5000000 -106.5000000 2759.5849609375 | cct -d 5 $PROJ_PIPELINE_7912

      38.50001      -106.50001    2744.47006           inf

Global DEMs

NASADEM