NASADEM - 3D Coordinate Reference System Transformation Resources

It’s reprocessed SRTM, so assume the same vertical reference as SRTM:

WGS (G1150) 3D + EGM1996 = EPSG:9055+5773

https://spatialreference.org/ref/epsg/9055/

https://spatialreference.org/ref/epsg/5773/

https://portal.opentopography.org/datasetMetadata?otCollectionID=OT.032021.4326.2

Also on Planetary Computer: https://planetarycomputer.microsoft.com/dataset/nasadem

%%bash

INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalinfo -nofl $INPUT

Driver: VRT/Virtual Raster
Size is 1288801, 421201
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            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",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-179.000138888888898,61.000138888888891)
Pixel Size = (0.000277777777778,-0.000277777777778)
Corner Coordinates:
Upper Left  (-179.0001389,  61.0001389) (179d 0' 0.50"W, 61d 0' 0.50"N)
Lower Left  (-179.0001389, -56.0001389) (179d 0' 0.50"W, 56d 0' 0.50"S)
Upper Right ( 179.0001389,  61.0001389) (179d 0' 0.50"E, 61d 0' 0.50"N)
Lower Right ( 179.0001389, -56.0001389) (179d 0' 0.50"E, 56d 0' 0.50"S)
Center      (   0.0000000,   2.5000000) (  0d 0' 0.00"E,  2d30' 0.00"N)
Band 1 Block=128x128 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

#%%bash
# Also inspect a single tile:
#GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
# gdalinfo /vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be/NASADEM_HGT_n38w107.tif

%%bash

INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -geoloc $INPUT -106.500 38.500

Report:
  Location: (261000P,81000L)
  Band 1:
    <LocationInfo><File>/vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be/NASADEM_HGT_n38w107.tif</File></LocationInfo>
    Value: 2762

Convert to Ellipsoid Height¶

Combined vertical shift grid + helmert transform OSGeo/PROJ#4362

%%bash

SSRS='EPSG:9055+5773'
TSRS='EPSG:7661'
PROJ_PIPELINE=`projinfo -s EPSG:9055+5773 -t EPSG:7661 -o PROJ --hide-ballpark -q --single-line`
INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt
OUTPUT=nasadem_7661.vrt

CPL_DEBUG=OFF PROJ_DEBUG=2 \
 PROJ_NETWORK=ON \
 GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalwarp -overwrite -wm 500 -r bilinear -ot Float32 -te -179.0001389 -56.0001389 179.0001389 61.0001389 -s_srs $SSRS -t_srs $TSRS -ct "${PROJ_PIPELINE}" ${INPUT} ${OUTPUT}

Creating output file that is 1288801P x 421201L.
Using internal nodata values (e.g. -32768) for image /vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt.
Copying nodata values from source /vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt to destination nasadem_7661.vrt.
Processing /vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

%%bash

INPUT=nasadem_7661.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -geoloc $INPUT -106.500 38.500

Report:
  Location: (261000P,81000L)
  Band 1:
    Value: 2747.1591796875

#projinfo -s EPSG:9055+5773 -t EPSG:7661 -o PROJ --hide-ballpark -q --single-line
#projinfo -s EPSG:9055 -t EPSG:7912 -o PROJ --hide-ballpark --single-line -q
pipeline='''+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1
  +step +proj=cart +ellps=WGS84
  +step +proj=helmert +x=-0.0007 +y=-0.0012 +z=0.0261 +rx=0 +ry=0 +rz=0
        +s=-0.00212 +dx=-0.0001 +dy=-0.0001 +dz=0.0019 +drx=0 +dry=0 +drz=0
        +ds=-0.00011 +t_epoch=2010 +convention=position_vector
  +step +inv +proj=cart +ellps=GRS80
  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
  +step +proj=axisswap +order=2,1'''
pipeline.replace('\n','')

'+proj=pipeline  +step +proj=axisswap +order=2,1  +step +proj=unitconvert +xy_in=deg +xy_out=rad  +step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1  +step +proj=cart +ellps=WGS84  +step +proj=helmert +x=-0.0007 +y=-0.0012 +z=0.0261 +rx=0 +ry=0 +rz=0        +s=-0.00212 +dx=-0.0001 +dy=-0.0001 +dz=0.0019 +drx=0 +dry=0 +drz=0        +ds=-0.00011 +t_epoch=2010 +convention=position_vector  +step +inv +proj=cart +ellps=GRS80  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m  +step +proj=axisswap +order=2,1'

%%bash

SSRS='EPSG:9055+5773'
TSRS='EPSG:7912'
PROJ_PIPELINE='+proj=pipeline  +step +proj=axisswap +order=2,1  +step +proj=unitconvert +xy_in=deg +xy_out=rad  +step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1  +step +proj=cart +ellps=WGS84  +step +proj=helmert +x=-0.0007 +y=-0.0012 +z=0.0261 +rx=0 +ry=0 +rz=0        +s=-0.00212 +dx=-0.0001 +dy=-0.0001 +dz=0.0019 +drx=0 +dry=0 +drz=0        +ds=-0.00011 +t_epoch=2010 +convention=position_vector  +step +inv +proj=cart +ellps=GRS80  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m  +step +proj=axisswap +order=2,1'
INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt
OUTPUT=nasadem_7912.vrt

CPL_DEBUG=OFF PROJ_DEBUG=2 \
 PROJ_NETWORK=ON \
 GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalwarp -overwrite -wm 500 -r bilinear -ot Float32 -te -179.0001389 -56.0001389 179.0001389 61.0001389 -s_srs $SSRS -t_srs $TSRS -ct "${PROJ_PIPELINE}" ${INPUT} ${OUTPUT}

Creating output file that is 1288801P x 421201L.
Using internal nodata values (e.g. -32768) for image /vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt.
Copying nodata values from source /vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt to destination nasadem_7912.vrt.
Processing /vsicurl/https://opentopography.s3.sdsc.edu/raster/NASADEM/NASADEM_be.vrt [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

%%bash

INPUT=nasadem_7912.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -geoloc $INPUT -106.500 38.500

Report:
  Location: (261000P,81000L)
  Band 1:
    Value: 2747.1611328125

Planetary Computer¶

Here we demonstrate using stac-geoparquet + GDAL GTI instead of a VRT

%%bash
# Feature Count: 14520, assets.elevation.href
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 VSICURL_PC_URL_SIGNING=YES \
 AZURE_STORAGE_ACCOUNT=pcstacitems \
 AZURE_STORAGE_SAS_TOKEN=`curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/pcstacitems/items | jq -r '.token'` \
 ogrinfo -al -so /vsiaz/items/nasadem.parquet

INFO: Open of `/vsiaz/items/nasadem.parquet'
      using driver `Parquet' successful.

Layer name: nasadem
Geometry: Polygon
Feature Count: 14520
Extent: (-179.000139, -56.000139) - (179.000139, 61.000139)
Layer SRS WKT:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    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]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Geometry Column = geometry
type: String (0.0)
stac_version: String (0.0)
stac_extensions: StringList (0.0)
id: String (0.0)
bbox: RealList (0.0)
links: String(JSON) (0.0)
assets.elevation.href: String (0.0)
assets.elevation.roles: StringList (0.0)
assets.elevation.title: String (0.0)
assets.elevation.type: String (0.0)
assets.rendered_preview.href: String (0.0)
assets.rendered_preview.rel: String (0.0)
assets.rendered_preview.roles: StringList (0.0)
assets.rendered_preview.title: String (0.0)
assets.rendered_preview.type: String (0.0)
assets.tilejson.href: String (0.0)
assets.tilejson.roles: StringList (0.0)
assets.tilejson.title: String (0.0)
assets.tilejson.type: String (0.0)
collection: String (0.0)
datetime: DateTime (UTC)
proj:bbox: RealList (0.0)
proj:epsg: Integer64 (0.0)
proj:shape: Integer64List (0.0)
proj:transform: RealList (0.0)

%%bash

CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 VSICURL_PC_URL_SIGNING=YES \
 AZURE_STORAGE_ACCOUNT=pcstacitems \
 AZURE_STORAGE_SAS_TOKEN=`curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/pcstacitems/items | jq -r '.token'` \
 gdalinfo GTI:/vsiaz/items/nasadem.parquet -oo LOCATION_FIELD=assets.elevation.href

Driver: GTI/GDAL Raster Tile Index
Files: none associated
Size is 1288802, 421202
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    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]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-179.000138889999988,61.000138890000002)
Pixel Size = (0.000277777777778,-0.000277777777778)
Corner Coordinates:
Upper Left  (-179.0001389,  61.0001389) (179d 0' 0.50"W, 61d 0' 0.50"N)
Lower Left  (-179.0001389, -56.0004167) (179d 0' 0.50"W, 56d 0' 1.50"S)
Upper Right ( 179.0004167,  61.0001389) (179d 0' 1.50"E, 61d 0' 0.50"N)
Lower Right ( 179.0004167, -56.0004167) (179d 0' 1.50"E, 56d 0' 1.50"S)
Center      (   0.0001389,   2.4998611) (  0d 0' 0.50"E,  2d29'59.50"N)
Band 1 Block=256x256 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

<VRTDataset rasterXSize="1296001" rasterYSize="417601">
  <SRS>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]</SRS>
  <GeoTransform> -1.8000013888888890e+02,  2.7777777777781459e-04,  0.0000000000000000e+00,  6.0000138888888891e+01,  0.0000000000000000e+00, -2.7777777777781459e-04</GeoTransform>
  <VRTRasterBand dataType="Int16" band="1">
    <NoDataValue>-32768</NoDataValue>

with open('nasadem.gti', 'w') as f:
    f.write("""<GDALTileIndexDataset>
    <IndexDataset>/vsiaz/items/nasadem.parquet</IndexDataset>
    <LocationField>assets.elevation.href</LocationField>
    <GeoTransform>-1.8000013888888890e+02,  2.7777777777781459e-04,0.0000000000000000e+00,6.0000138888888891e+01,0.0000000000000000e+00,-2.7777777777781459e-04</GeoTransform>
    <XSize>1296001</XSize>
    <YSize>417601</YSize>
    <BandCount>1</BandCount>
    <DataType>Int16</DataType>
</GDALTileIndexDataset>
""")

%%bash

CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 VSICURL_PC_URL_SIGNING=YES \
 AZURE_STORAGE_ACCOUNT=pcstacitems \
 AZURE_STORAGE_SAS_TOKEN=`curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/pcstacitems/items | jq -r '.token'` \
 gdalinfo nasadem.gti

Driver: GTI/GDAL Raster Tile Index
Files: nasadem.gti
Size is 1296001, 417601
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    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]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000138888888898,60.000138888888891)
Pixel Size = (0.000277777777778,-0.000277777777778)
Corner Coordinates:
Upper Left  (-180.0001389,  60.0001389) (180d 0' 0.50"W, 60d 0' 0.50"N)
Lower Left  (-180.0001389, -56.0001389) (180d 0' 0.50"W, 56d 0' 0.50"S)
Upper Right ( 180.0001389,  60.0001389) (180d 0' 0.50"E, 60d 0' 0.50"N)
Lower Right ( 180.0001389, -56.0001389) (180d 0' 0.50"E, 56d 0' 0.50"S)
Center      (   0.0000000,   2.0000000) (  0d 0' 0.00"E,  2d 0' 0.00"N)
Band 1 Block=256x256 Type=Int16, ColorInterp=Undefined

%%bash

INPUT=nasadem.gti
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 VSICURL_PC_URL_SIGNING=YES \
 AZURE_STORAGE_ACCOUNT=pcstacitems \
 AZURE_STORAGE_SAS_TOKEN=`curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/pcstacitems/items | jq -r '.token'` \
 gdallocationinfo -geoloc $INPUT -106.500 38.500

Report:
  Location: (264600P,77400L)
  Band 1:
    <LocationInfo><File>/vsicurl/https://nasademeuwest.blob.core.windows.net/nasadem-cog/v001/NASADEM_HGT_n38w107.tif</File></LocationInfo>
    Value: 2762

%%bash
# Query Tif directly
INPUT=/vsicurl/https://nasademeuwest.blob.core.windows.net/nasadem-cog/v001/NASADEM_HGT_n38w107.tif
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 VSICURL_PC_URL_SIGNING=YES \
 gdallocationinfo -geoloc $INPUT -106.500 38.500

Report:
  Location: (1800P,1800L)
  Band 1:
    Value: 2762


%%bash

SSRS='EPSG:9055+5773'
TSRS='EPSG:7912'
PROJ_PIPELINE='+proj=pipeline  +step +proj=axisswap +order=2,1  +step +proj=unitconvert +xy_in=deg +xy_out=rad  +step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1  +step +proj=cart +ellps=WGS84  +step +proj=helmert +x=-0.0007 +y=-0.0012 +z=0.0261 +rx=0 +ry=0 +rz=0        +s=-0.00212 +dx=-0.0001 +dy=-0.0001 +dz=0.0019 +drx=0 +dry=0 +drz=0        +ds=-0.00011 +t_epoch=2010 +convention=position_vector  +step +inv +proj=cart +ellps=GRS80  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m  +step +proj=axisswap +order=2,1'
INPUT=nasadem.gti
OUTPUT=nasadem_7912.gti.vrt

CPL_DEBUG=OFF PROJ_DEBUG=2 \
 PROJ_NETWORK=ON \
 GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 VSICURL_PC_URL_SIGNING=YES \
 AZURE_STORAGE_ACCOUNT=pcstacitems \
 AZURE_STORAGE_SAS_TOKEN=`curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/pcstacitems/items | jq -r '.token'` \
 gdalwarp -overwrite -wm 500 -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 -r bilinear -ot Float32 -te -180.0001389 -56.0001389 180.0001389 60.0001389 -s_srs $SSRS -t_srs $TSRS -ct "${PROJ_PIPELINE}" ${INPUT} ${OUTPUT}

ERROR 1: Point outside of projection domain

Creating output file that is 1296001P x 417601L.
Processing nasadem.gti [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

%%bash
# Does this match opentopography VRT? Value: 2747.1611328125, yes
INPUT=nasadem_7912.gti.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 VSICURL_PC_URL_SIGNING=YES \
 AZURE_STORAGE_ACCOUNT=pcstacitems \
 AZURE_STORAGE_SAS_TOKEN=`curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/pcstacitems/items | jq -r '.token'` \
 gdallocationinfo -geoloc $INPUT -106.500 38.500

Report:
  Location: (264600P,77400L)
  Band 1:
    Value: 2747.1611328125

%%bash
# Does this match using PROJ directly?
INPUT=nasadem_7912.gti.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 VSICURL_PC_URL_SIGNING=YES \
 AZURE_STORAGE_ACCOUNT=pcstacitems \
 AZURE_STORAGE_SAS_TOKEN=`curl -s https://planetarycomputer.microsoft.com/api/sas/v1/token/pcstacitems/items | jq -r '.token'` \
 gdallocationinfo -r bilinear -geoloc $INPUT -106.500000 38.500000

Report:
  Location: (264600.500023667P,77400.5000251724L)
  Band 1:
    Value: 2747.16103365231

%%bash
# 7912 Single point using proj, agrees w/ GDAL to w/n mm
PROJ_PIPELINE='+proj=pipeline  +step +proj=axisswap +order=2,1  +step +proj=unitconvert +xy_in=deg +xy_out=rad  +step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1  +step +proj=cart +ellps=WGS84  +step +proj=helmert +x=-0.0007 +y=-0.0012 +z=0.0261 +rx=0 +ry=0 +rz=0        +s=-0.00212 +dx=-0.0001 +dy=-0.0001 +dz=0.0019 +drx=0 +dry=0 +drz=0        +ds=-0.00011 +t_epoch=2010 +convention=position_vector  +step +inv +proj=cart +ellps=GRS80  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m  +step +proj=axisswap +order=2,1'
#USE proj directly to convert this point
echo 38.5000000 -106.5000000 2762.0 | cct -d 5 $PROJ_PIPELINE

      38.50000      -106.50000    2747.16309           inf