Copernicus DEM - 3D Coordinate Reference System Transformation Resources

This notebook illustrates a way to encode a 3D proj transformation into a GDAL VRT so that you can load elevation ellipsoid heights rather than geoid heights. This can happen “on-the-fly”, so there is no need to download and save local files.

We’ll use: https://portal.opentopography.org/raster?opentopoID=OTSDEM.032021.4326.3

Horizontal Coordinates: WGS84 [EPSG: 4326]
Vertical Coordinates: EGM2008 [EPSG: 3855]

In GDAL we can specify a full 3D CRS as a compound CRS EPSG:9055+3855

%%bash
# NOTE: 'GEOCRS' SCOPE["Horizontal component of 3D system."], CS[ellipsoidal,2],
# ENSEMBLE["World Geodetic System 1984 ensemble",
# ENSEMBLEACCURACY[2.0]],
gdalsrsinfo EPSG:4326


PROJ.4 : +proj=longlat +datum=WGS84 +no_defs

OGC WKT2:2019 :
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]]

%%bash
# NOTE: GEOGCRS["WGS 84 (G1150)" + VERTCRS["EGM2008 height"
projinfo EPSG:9055+3855

PROJ.4 string:

Cannot open https://cdn.proj.org/Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz: HTTP error 404: <?xml version="1.0" encoding="UTF-8"?>
<Error><Code>NoSuchKey</Code><Message>The specified key does not exist.</Message><Key>Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz</Key><RequestId>G5EQJJJ1G3SNACPW</RequestId><HostId>3mKkT5qXDZdI3B0lNvB9btUiJPPU200jcVJOcrk5Vxsf4iXklsUdHMyJwrQ3l1g4UswitgLAgvU=</HostId></Error>
Cannot open https://cdn.proj.org/Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz: HTTP error 404: <?xml version="1.0" encoding="UTF-8"?>
<Error><Code>NoSuchKey</Code><Message>The specified key does not exist.</Message><Key>Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz</Key><RequestId>G5EQJJJ1G3SNACPW</RequestId><HostId>3mKkT5qXDZdI3B0lNvB9btUiJPPU200jcVJOcrk5Vxsf4iXklsUdHMyJwrQ3l1g4UswitgLAgvU=</HostId></Error>

+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +geoidgrids=us_nga_egm08_25.tif +geoid_crs=WGS84 +vunits=m +no_defs +type=crs

WKT2:2019 string:
COMPOUNDCRS["WGS 84 (G1150) + EGM2008 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["EGM2008 height",
        VDATUM["EGM2008 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",3855]]]

Convert to ellipsoid height¶

We are interested in comparing DEM elevations to altimeter elevations made by ICESAT-2. Which uses EPSG:7912 (ITRF2014 ellipsoid heights).

The simplest transform to ellipsoid height would be to keep the same datum realization and apply a vertical shift to account for the geoid offset:

Caution: use 3D CRS¶

%%bash

# CAREFUL! 3D -> 2D definitions just ignore the Z-axis by default.
projinfo -s EPSG:9055+3855 -t EPSG:9055 -o PROJ --hide-ballpark

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

unknown id, Null geographic offset from WGS 84 (G1150) to WGS 84 (G1150), 0 m, World.

PROJ string:
+proj=noop

%%bash

# WGS (G1150) 3D + EGM2008 -> WGS (G1150) 3D
# NOTE: +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
# NOTE: The shift grid metadata uses WGS 84 ENSEMBLE (EPSG:4979) to EGM2008 height (EPSG:3855), resulting 5m accuracy?
projinfo -s EPSG:9055+3855 -t EPSG:7661 -o PROJ --hide-ballpark

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

unknown id, Inverse of WGS 84 to EGM2008 height (1) using WGS 84 to WGS 84 (G1150), 5 m, World.

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

%%bash
# To get a properly formatted single-line string to pass to other programs like gdalwarp:
projinfo -s EPSG:9055+3855 -t EPSG:7661 -o PROJ --hide-ballpark --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_nga_egm08_25.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

# If you're going to use this grid a lot, download it!
!projsync --file us_nga_egm08_25

Downloading from https://cdn.proj.org into /Users/scotthenderson/Library/Application Support/proj
https://cdn.proj.org/us_nga_egm08_25.tif already downloaded.

!gdalinfo /vsicurl/https://cdn.proj.org/us_nga_egm08_25.tif

Driver: GTiff/GeoTIFF
Files: /vsicurl/https://cdn.proj.org/us_nga_egm08_25.tif
       /vsicurl/https://cdn.proj.org/us_nga_egm08_25.tif.aux.xml
Size is 8640, 4321
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.020833333333343,90.020833333333329)
Pixel Size = (0.041666666666667,-0.041666666666667)
Metadata:
  TIFFTAG_IMAGEDESCRIPTION=WGS 84 (EPSG:4979) to EGM2008 height (EPSG:3855). Converted from egm08_25.gtx (last modified at 2018/10/08)
  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=3855
  TYPE=VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL
  AREA_OR_POINT=Point
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
  PREDICTOR=3
Corner Coordinates:
Upper Left  (-180.0208333,  90.0208333) (180d 1'15.00"W, 90d 1'15.00"N)
Lower Left  (-180.0208333, -90.0208333) (180d 1'15.00"W, 90d 1'15.00"S)
Upper Right ( 179.9791667,  90.0208333) (179d58'45.00"E, 90d 1'15.00"N)
Lower Right ( 179.9791667, -90.0208333) (179d58'45.00"E, 90d 1'15.00"S)
Center      (  -0.0208333,   0.0000000) (  0d 1'15.00"W,  0d 0' 0.00"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = geoid_undulation
  Min=-50.018 Max=85.824 
  Unit Type: metre
  Metadata:
    STATISTICS_APPROXIMATE=YES
    STATISTICS_MAXIMUM=85.823783874512
    STATISTICS_MEAN=1.0811106698347
    STATISTICS_MINIMUM=-50.018199920654
    STATISTICS_STDDEV=26.622038806097
    STATISTICS_VALID_PERCENT=100

Helmert transform to more recent WGS/ITRF Realization¶

The helmert transform accounts for updates to the datum realization and plate motion over time (so t_epoch is specified). This provides a matrix-transform to go from one WGS84/ITRF realization to another. These transforms are applied to ECEF coordinates. You can get parameters via online calculators (https://itrf.ign.fr/en/solutions/transformations) or look them up in the PROJ database:

%%bash
# No grid shift but next WGS realization
# WGS (G1150) 3D + EGM2008 -> WGS (G1674) 3D
projinfo -s EPSG:7661 -t EPSG:7663 -o PROJ

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

unknown id, Conversion from WGS 84 (G1150) (geog3D) to WGS 84 (G1150) (geocentric) + WGS 84 (G1150) to WGS 84 (G1674) (1) + Conversion from WGS 84 (G1674) (geocentric) to WGS 84 (G1674) (geog3D), 0.02 m, World

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
  +step +proj=cart +ellps=WGS84
  +step +proj=helmert +x=-0.0024 +y=0.0016 +z=0.0232 +rx=-0.00027 +ry=0.00027
        +rz=-0.00038 +s=0.00208 +dx=-0.0001 +dy=-0.0001 +dz=0.0018 +drx=0 +dry=0
        +drz=0 +ds=-8e-05 +t_epoch=2005 +convention=coordinate_frame
  +step +inv +proj=cart +ellps=WGS84
  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
  +step +proj=axisswap +order=2,1

!projinfo -s EPSG:7661 -t EPSG:7912 -o PROJ --single-line -q

+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +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
# Combined NOTE: this requires PROJ>=9.6
# https://github.com/OSGeo/PROJ/issues/4362
#projinfo -s EPSG:9055+3855 -t EPSG:7663 -o proj

    +proj=pipeline
      +step +proj=axisswap +order=2,1
      +step +proj=unitconvert +xy_in=deg +xy_out=rad
      +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
      +step +proj=cart +ellps=WGS84
      +step +proj=helmert +x=-0.0024 +y=0.0016 +z=0.0232 +rx=-0.00027 +ry=0.00027
            +rz=-0.00038 +s=0.00208 +dx=-0.0001 +dy=-0.0001 +dz=0.0018 +drx=0 +dry=0
            +drz=0 +ds=-8e-05 +t_epoch=2005 +convention=coordinate_frame
      +step +inv +proj=cart +ellps=WGS84
      +step +proj=unitconvert +xy_in=rad +xy_out=deg
      +step +proj=axisswap +order=2,1

# Note: using '--spatial-test intersects' would bring more results (62)
# See https://proj.org/en/stable/apps/projinfo.html#cmdoption-projinfo-spatial-test
# For global grid shifts --spatial-test interests isn't ideal because PROJ will offer options of using regional datums like GDA2020 (defined just for Australia)
#!projinfo -s EPSG:9055+3855 -t EPSG:7912 -o proj

+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_nga_egm08_25.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 +xy_out=deg
  +step +proj=axisswap +order=2,1

Encoding transform in a GDAL VRT¶

%%bash

INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt

# NOTE: -nofl is for "no file list" (Useful For VRTs with tons of rasters!)
# NOTE: EPSG:4326 is wrong :( so we need to override it!

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalinfo -nofl $INPUT

Driver: VRT/Virtual Raster
Size is 1296001, 626401
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.000138899999996,84.000138899999996)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  (-180.0001389,  84.0001389) (180d 0' 0.50"W, 84d 0' 0.50"N)
Lower Left  (-180.0001389, -90.0001389) (180d 0' 0.50"W, 90d 0' 0.50"S)
Upper Right ( 180.0001389,  84.0001389) (180d 0' 0.50"E, 84d 0' 0.50"N)
Lower Right ( 180.0001389, -90.0001389) (180d 0' 0.50"E, 90d 0' 0.50"S)
Center      (  -0.0000000,  -3.0000000) (  0d 0' 0.00"W,  3d 0' 0.00"S)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray

%%bash

# A way to trim out of bounds pixels (< -90 latitude) in theory this will get rid of
# ERROR 1: PROJ: vgridshift: Invalid latitude
# But would need to adjust transform on the VRT as well...

GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalinfo -nofl "vrt:///tmp/COP30_hh.vrt?srcwin=0,-2,1296001,626399"

Driver: VRT/Virtual Raster
Size is 1296001, 626399
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.000138899999996,84.000694455555546)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  (-180.0001389,  84.0006945) (180d 0' 0.50"W, 84d 0' 2.50"N)
Lower Left  (-180.0001389, -89.9990278) (180d 0' 0.50"W, 89d59'56.50"S)
Upper Right ( 180.0001389,  84.0006945) (180d 0' 0.50"E, 84d 0' 2.50"N)
Lower Right ( 180.0001389, -89.9990278) (180d 0' 0.50"E, 89d59'56.50"S)
Center      (  -0.0000000,  -2.9991667) (  0d 0' 0.00"W,  2d59'57.00"S)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray

%%bash

# ENSURE EGM2008 applied

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_egm08_25.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1'
echo $PROJ_PIPELINE
INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt
OUTPUT=COP30_hh_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 \
    -te -180.000138888888898 -90.000138888888891 180.000138888888891 84.000138888888891 \
    -tr 0.000277777777778 0.000277777777778 \
    -s_srs EPSG:9055+3855 -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_egm08_25.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 626401L.
Processing /vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

%%bash

# Use GDAL Warp to *encode* the appropriate PROJ pipeline
# Pass extents explicity using PROJ pipeline
# https://github.com/OSGeo/gdal/issues/11610
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_egm08_25.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 +xy_out=deg +step +proj=axisswap +order=2,1'
echo $PROJ_PIPELINE
INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt
OUTPUT=/tmp/COP30_hh_7912.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 \
    -te -180.000138888888898 -90.000138888888891 180.000138888888891 84.000138888888891 \
    -tr 0.000277777777778 0.000277777777778 \
    -s_srs EPSG:9055+3855 -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_egm08_25.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 +xy_out=deg +step +proj=axisswap +order=2,1
Creating output file that is 1296001P x 626401L.
Processing /vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

Validation¶

Query points & compare to transforming single points via PROJ
Transform a single tiff from within the VRT and compare to on-the-fly transform of same subset region

%%bash

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

Report:
  Location: (264600.50004P,163800.50004L)
  Band 1:
    <LocationInfo><File>/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh/Copernicus_DSM_10_N38_00_W107_00_DEM.tif</File></LocationInfo>
    Value: 2758.54790051775

%%bash
# Sample exactly the same point from the embedded TIFF
# NOTE: not sure why these don't agree at the sub-mm level...
INPUT=/vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh/Copernicus_DSM_10_N38_00_W107_00_DEM.tif
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -r bilinear -geoloc $INPUT -106.500 38.500

Report:
  Location: (1800.5P,1800.5L)
  Band 1:
    Value: 2758.5478515625

%%bash
# Confirm Geoid shift ~10m yep.
INPUT=/tmp/COP30_hh_7912.vrt
CPL_DEBUG=OFF GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdallocationinfo -r bilinear -geoloc $INPUT -106.500 38.500

ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude
ERROR 1: PROJ: vgridshift: Invalid latitude

Report:
  Location: (264600.499999788P,163800.499999869L)
  Band 1:
    Value: 2744.04248051403

%%bash
# 7912 Single point using proj, agrees w/ GDAL to sub-mm levels
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_egm08_25.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 +xy_out=deg +step +proj=axisswap +order=2,1'
#USE proj directly to convert this point
echo 38.5000000 -106.5000000 2758.54790051775 | cct -d 5 $PROJ_PIPELINE

      38.50000      -106.50000    2744.04248           inf