[gdal-dev] gdalwarp vertical datum transformation with +geoidgrids with no changes

[Karl Zinglersen]

I try to perform a vertical transformation from WGS84 ellipsoid heights to
geoid MSL heights on a Geotiff with a geoid GTX raster translated from
Geotiff.

gdalwarp -overwrite -s_srs "+proj=longlat +datum=WGS84 +no_defs
+geoidgrids=C:\TEMP\gr20006g06.gtx" -t_srs EPSG:4326 -r bilinear -of GTiff
C:/TEMP/gimdem1_1_epsg4326.tif C:/TEMP/gimdem1_1_epsg4326_geoid.tif

The transformation runs, but there are no changes in cell values of the
output raster as expected.
I can do a gdal_calc between the input raster and the geoid raster to an
output raster, which changes the pixel values, but requires a number of
previous steps: align, clip, warp, resampling of the geoid raster.

I would think the +geoidgrids warping is the most correct?

How can I change my data or my gdalwarp flags and inputs to do a
transformation which actually use the geoid grid?

Some information on the input DEM and geoid:

Input DEM:
Driver: GTiff/GeoTIFF
Files: C:/TEMP/gimdem1_1_epsg4346.tif
Size is 14727, 9815
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-54.040095691152260,67.573648952718571)
Pixel Size = (0.000424686603675,-0.000424686603675)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -54.0400957,  67.5736490) ( 54d 2'24.34"W, 67d34'25.14"N)
Lower Left  ( -54.0400957,  63.4053499) ( 54d 2'24.34"W, 63d24'19.26"N)
Upper Right ( -47.7857361,  67.5736490) ( 47d47' 8.65"W, 67d34'25.14"N)
Lower Right ( -47.7857361,  63.4053499) ( 47d47' 8.65"W, 63d24'19.26"N)
Center      ( -50.9129159,  65.4894994) ( 50d54'46.50"W, 65d29'22.20"N)
Band 1 Block=14727x1 Type=Int16, ColorInterp=Gray


Geoid raster:
Driver: GTX/NOAA Vertical Datum .GTX
Files: C:/TEMP/gr20006g06.gtx
       C:/TEMP/gr20006g06.gtx.aux.xml
Size is 651, 251
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9108"]],
    AUTHORITY["EPSG","4326"]]
Origin = (-75.000000000000000,84.000000000000000)
Pixel Size = (0.100000000000000,-0.100000000000000)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  ( -75.0000000,  84.0000000) ( 75d 0' 0.00"W, 84d 0' 0.00"N)
Lower Left  ( -75.0000000,  58.9000000) ( 75d 0' 0.00"W, 58d54' 0.00"N)
Upper Right (  -9.9000000,  84.0000000) (  9d54' 0.00"W, 84d 0' 0.00"N)
Lower Right (  -9.9000000,  58.9000000) (  9d54' 0.00"W, 58d54' 0.00"N)
Center      ( -42.4500000,  71.4500000) ( 42d27' 0.00"W, 71d27' 0.00"N)
Band 1 Block=651x1 Type=Float32, ColorInterp=Undefined
  NoData Value=9999






-- 
Karl Brix Zinglersen
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20160308/01169e09/attachment.html>