3DEP Seamless DEMs

%%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

%%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

# 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

# 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')

# 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

# 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

# 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

%%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}

%%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

%%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