Skip to article frontmatterSkip to article content

SRTM 30m

Survey Date: 02/11/2000 - 02/22/2000
Horizontal Coordinates: WGS84 [EPSG: 4326]
Vertical Coordinates: WGS84 (EGM96 GEOID)

https://portal.opentopography.org/raster?opentopoID=OTSRTM.082015.4326.1

NOTE: 4326 is not ideal, instead use specific Geoid reference:

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

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

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

Warning: the opentopo SRTM data is also not very cloud optimized see https://lists.osgeo.org/pipermail/gdal-dev/2025-January/thread.html

%%bash

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalinfo -nofl /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt
Driver: VRT/Virtual Raster
Size is 1296001, 417601
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 = (-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=128x128 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

VRT Connection Strings

a convenient way to override the CRS in metadata is with vrt connection strings:

https://gdal.org/en/stable/drivers/raster/vrt.html#vrt-connection-string

%%bash

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalinfo -nofl vrt:///vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt?a_srs=EPSG:9055+5773
Driver: VRT/Virtual Raster
Size is 1296001, 417601
Coordinate System is:
COMPOUNDCRS["WGS 84 (G1150) + EGM96 height",
    GEOGCRS["WGS 84 (G1150)",
        DYNAMIC[
            FRAMEEPOCH[2001]],
        DATUM["World Geodetic System 1984 (G1150)",
            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]],
        USAGE[
            SCOPE["Geodesy. Navigation and positioning using GPS satellite system."],
            AREA["World."],
            BBOX[-90,-180,90,180]],
        ID["EPSG",9055]],
    VERTCRS["EGM96 height",
        VDATUM["EGM96 geoid"],
        CS[vertical,1],
            AXIS["gravity-related height (H)",up,
                LENGTHUNIT["metre",1]],
        USAGE[
            SCOPE["Geodesy."],
            AREA["World."],
            BBOX[-90,-180,90,180]],
        ID["EPSG",5773]]]
Data axis to CRS axis mapping: 2,1,3
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=128x128 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

%%bash

# NOTE: also SRTM_GL1_Ellip/SRTM_GL1_Ellip_srtm.vrt !
# which looks like geoid -> ellipsoid precomputed per tile (SRTM_GL1_Ellip_srtm/North/North_30_60/N52E119_wgs84.tif)
GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalinfo -nofl /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1_Ellip/SRTM_GL1_Ellip_srtm.vrt
Driver: VRT/Virtual Raster
Size is 1296001, 417601
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 = (-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=128x128 Type=Float32, ColorInterp=Gray
  NoData Value=-32768

%%bash
# Look at a single TIF w/n the VRT

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalinfo /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm/N38W107.tif
Driver: GTiff/GeoTIFF
Files: /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm/N38W107.tif
Size is 3601, 3601
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 = (-107.000138888888884,39.000138888888891)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  OVR_RESAMPLING_ALG=NEAREST
  AREA_OR_POINT=Point
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-107.0001389,  39.0001389) (107d 0' 0.50"W, 39d 0' 0.50"N)
Lower Left  (-107.0001389,  37.9998611) (107d 0' 0.50"W, 37d59'59.50"N)
Upper Right (-105.9998611,  39.0001389) (105d59'59.50"W, 39d 0' 0.50"N)
Lower Right (-105.9998611,  37.9998611) (105d59'59.50"W, 37d59'59.50"N)
Center      (-106.5000000,  38.5000000) (106d30' 0.00"W, 38d30' 0.00"N)
Band 1 Block=512x512 Type=Int16, ColorInterp=Gray
  NoData Value=-32768
  Overviews: 1801x1801, 901x901, 451x451

Projinfo to check possible CRS Transforms

# Convert to ellipsoid height
!projinfo -s EPSG:9055+5773 -t EPSG:7661 -o PROJ --hide-ballpark -q --single-line
+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=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

# Convert to ITRF2014  (ICESat-2)
!projinfo -s EPSG:9055 -t EPSG:7912 -o PROJ --hide-ballpark
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Conversion from WGS 84 (G1150) (geog2D) to WGS 84 (G1150) (geocentric) + WGS 84 (G1150) to ITRF2000 (1) + ITRF2000 to ITRF2014 (1) + Conversion from ITRF2014 (geocentric) to ITRF2014 (geog3D), 0.03 m, World

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +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

# Combined (inject shift grid https://github.com/OSGeo/PROJ/issues/4362
# Pipeline as a single line to pass to GDAL
#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'
!GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR gdalinfo /vsicurl/https://cdn.proj.org/us_nga_egm96_15.tif
Driver: GTiff/GeoTIFF
Files: /vsicurl/https://cdn.proj.org/us_nga_egm96_15.tif
Size is 1440, 721
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,3],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Geodesy. Navigation and positioning using GPS satellite system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4979]]
Data axis to CRS axis mapping: 2,1,3
Origin = (-180.125000000000000,90.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Metadata:
  TIFFTAG_IMAGEDESCRIPTION=WGS 84 (EPSG:4979) to EGM96 height (EPSG:5773). Converted from egm96_15.gtx (last modified at 2018/10/05)
  TIFFTAG_DATETIME=2019:12:27 00:00:00
  TIFFTAG_COPYRIGHT=Derived from work by NGA. Public Domain
  area_of_use=World
  target_crs_epsg_code=5773
  TYPE=VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL
  AREA_OR_POINT=Point
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
  PREDICTOR=3
Corner Coordinates:
Upper Left  (-180.1250000,  90.1250000) (180d 7'30.00"W, 90d 7'30.00"N)
Lower Left  (-180.1250000, -90.1250000) (180d 7'30.00"W, 90d 7'30.00"S)
Upper Right ( 179.8750000,  90.1250000) (179d52'30.00"E, 90d 7'30.00"N)
Lower Right ( 179.8750000, -90.1250000) (179d52'30.00"E, 90d 7'30.00"S)
Center      (  -0.1250000,   0.0000000) (  0d 7'30.00"W,  0d 0' 0.01"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = geoid_undulation
  Unit Type: metre

Encode Specific PROJ transform in a VRT

%%bash
# JUST vertical shift grid
# Copied combined pipeline from above
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=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1'
echo $PROJ_PIPELINE

# NOTE: enforce bounds and size to match original global VRT geotransform
# https://github.com/OSGeo/gdal/issues/11610

# Size is 1296001, 417601
#Pixel Size = (0.000277777777778,-0.000277777777778)
#Origin = (-180.000138888888898,60.000138888888891)
#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)
# -te -180.000138888888898 -56.000138888888891 180.000138888888891 60.000138888888891
# -tr 0.000277777777778 0.000277777777778
# -ts 1296001 417601
# -of Float32 (otherwise keeps Int16 of original data)
# -r nearest (for sampling original TIFs since the grid is the same, but sampling us_nga_egm96_15.tif uses bilinear interpolation)

#NOTE: at least for point queries, setting minmal VRT blocksize of 128x128 seems to result in fewest GET requests

INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt
OUTPUT=/tmp/SRTM_GL1_srtm_7661.vrt

CPL_DEBUG=OFF PROJ_DEBUG=2 \
 PROJ_NETWORK=ON \
 GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalwarp -overwrite -wm 500 -co BLOCKXSIZE=128 -co BLOCKYSIZE=128 \
    -ot Float32 \
    -te -180.000138888888898 -56.000138888888891 180.000138888888891 60.000138888888891 \
    -ts 1296001 417601 \
    -r nearest \
    -s_srs EPSG:9055+5773 -t_srs EPSG:7661 -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_nga_egm96_15.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
Creating output file that is 1296001P x 417601L.
Using internal nodata values (e.g. -32768) for image /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt.
Copying nodata values from source /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt to destination /tmp/SRTM_GL1_srtm_7661.vrt.
Processing /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

%%bash
# Vertical shift + Helmert
# Copied combined pipeline from above
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'
echo $PROJ_PIPELINE

INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt
OUTPUT=/tmp/SRTM_GL1_srtm_7912.vrt

CPL_DEBUG=OFF PROJ_DEBUG=3 \
 PROJ_NETWORK=ON \
 GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalwarp -wm 500 -overwrite \
    -co BLOCKXSIZE=128 -co BLOCKYSIZE=128 \
    -ot Float32 \
    -te -180.000138888888898 -56.000138888888891 180.000138888888891 60.000138888888891 \
    -ts 1296001 417601 \
    -r nearest \
    -s_srs EPSG:9055+5773 -t_srs EPSG:7912 -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_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
Creating output file that is 1296001P x 417601L.
Using internal nodata values (e.g. -32768) for image /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt.
Copying nodata values from source /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt to destination /tmp/SRTM_GL1_srtm_7912.vrt.
Processing /vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

Verify that VRT gives expected values

Because Opentopography has pre-computed tifs with the geoid -> ellipsoid transform, we can confirm that those have the same elevation values as using the VRT we just created (which computes the transform on-the-fly).

7661 and 7912 should agree w/n centimeters

%%bash

# Point query (nearest / point-containing pixel value) GEOID HEIGHT

INPUT=/tmp/SRTM_GL1_srtm_7912.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -geoloc $INPUT -106.500 38.500
Report:
  Location: (264600P,77400L)
  Band 1:
    Value: 2751.1630859375

%%bash

# Point query (nearest / point-containing pixel value) GEOID HEIGHT
# Hmmm, or maybe those ranges reflect the coverage of 5 EGM96 pixels to do bilinear interpolation?
INPUT=/tmp/SRTM_GL1_srtm_7661.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -geoloc $INPUT -106.500 38.500
Report:
  Location: (264600P,77400L)
  Band 1:
    Value: 2751.1591796875

%%bash

# Check against pre-computed value... not exactly sure how they generate these files, but should be close to above values

INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1_Ellip/SRTM_GL1_Ellip_srtm.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -geoloc $INPUT -106.500 38.500
Report:
  Location: (264600P,77400L)
  Band 1:
    <LocationInfo><File>/vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1_Ellip/SRTM_GL1_Ellip_srtm/North/North_30_60/N38W107_wgs84.tif</File></LocationInfo>
    Value: 2751.1591796875

%%bash

# Verify geoid height is offset from Ellipsoid Height

INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -geoloc $INPUT -106.500 38.500
Report:
  Location: (264600P,77400L)
  Band 1:
    <LocationInfo><File>/vsicurl/https://opentopography.s3.sdsc.edu/raster/SRTM_GL1/SRTM_GL1_srtm/N38W107.tif</File></LocationInfo>
    Value: 2766

Taking it further:

  • Verify this is working as expected by looking at validation points around the globe
Global DEMs
Copernicus DEM
Global DEMs
NASADEM