gtx file header has invalid extents

[jolanglois]

Hello,
I have converted NRCAN's grid shift files to gtx files using your tool byn2gtx. When I tried to use those gtx files to convert a .tif from EPSG:4326 to EPSG4326+5713 with gdal (installed on Windows 7 with Osgeo4W) I had this error message :

gtx file header has invalid extents, corrupt?
ERROR 1: failed to load datum shift file
ERROR 1: GDALSuggestedWarpOutput() failed because the passed tranformer failed.

I took the byn files on NRCan web site: http://webapp.geod.nrcan.gc.ca/geod/data-donnees/geoid.php?locale=fr

Here is the message I had when using byn2gtx (seems like there is effectively a problem with the extent in the gtx.):

Could you help me with this issue?
Thanks!

[rskelly]

Hi, it looks like the spec for byn files has changed since I wrote this. Let me have a look.

[jolanglois]

Ok thanks!
I had a look at the code and when I remove the double_flip function before printing the gtx header (in gtx.c, gtx_print_header function) I have good values for lower-left Corner and Deltas in the console, but I still have wrong values written in the gtx file.

[rskelly]

Yeah, I noticed that too. There are a couple of new fields in the header too.

[jolanglois]

Do you still have the gtx files you converted with this tool (CGG2013n83.gtx and HT2_0.gtx) If so, could you made them available?

[rskelly]

http://dijital.ca/files/nrcan_geoids.zip

These were produced a couple of years ago, but they work fine in proj4, etc.

[jolanglois]

Thanks for the gtx files!! Now I'm trying to use them to convert .las with las2las and .tif with gdalwarp and it failed.

gdalwarp -s_srs EPSG:4617+5713 -t_srs EPSG:4617 D:\Inuvik_CGVD28.tif D:\Inuvik_ellipsoid.tif
Output tif elevation stays exactly the same as the input even if this zone is covered by the gtx file.

las2las -i C:\Inuvik\CGVD28.las -o C:\Inuvik\CGVD28_set.las --a_srs '+proj=latlong +ellps=GRS80 +geoidgrids=HT2_0.gtx' --t_srs EPSG:4326
I got the error message : error: multiple occurrences.

I know it might not be the best place to ask for this, but I can't find any documentation on this. Do you think you could help me with those issues?

[rskelly]

To use the compound SRID you have to add it to your vertcs.override.csv file, which is in GDAL's share folder (I'm not sure where it is on Windows.) There're already records for 6647 and 5713 in vertcs.csv, but they don't include the geoidgrids.

las2las works with the compound CRS but gdalwarp doesn't. It won't apply geoid shifts, as far as I know.

I'm not familiar with the error you got, but you could try: "--a_srs epsg:4326+5713 --t_srs epsg:4326", since it looks like you're just trying to back out the geoid undulation to get ellipsoidal heights...

[rskelly]

I've committed fixes that should work with the new BYN spec. Also respects the correct nodata value for GTX files.

[jolanglois]

Thanks! It works fine. I also manage to convert vertical datum with las2las using EPSG:4326. I am just wondering why the gtx file Datum is set to WGS84 while the byn file Datum is set to NAD83 (ellipsoid GRS80). I know the difference between the two ellipsoid is millimetric so it is not a big issue, but I would like to know why you're not writting the original datum (GRS80) in the gtx file?

[rskelly]

I guess it's because Canada's official datum is NAD83, and it's fixed to the North American plate, which rotates over time. WGS84 is a geocentric reference frame, equivalent to ITRF. When using NAD83(CSRS) you have the concept of the epoch, and coordinates change over time relative to the reference frame. I noticed there is now space in the byn header for the epoch and reference frame.

Also, I was going to mention that if you want to use the geoid model with gdalwarp, you can resample it to the resolution of the raster you want to shift, and then just add or subtract them with gdal_calc.py.

[jolanglois]

Does it mean that there is no tool available in gdal to do vertical datum conversion of a raster DEM. I'm looking for a single command line to execute the vertical datum conversion of a raster like the las2las of Liblas does for Point Cloud data.

[rskelly]

I don't think so. That would be a good question for gis.stackexchange.com. My guess would be that there's no single tool because it's outside of the mandate of gdalwarp, and can easily be done by gdal_calc.

[jolanglois]

I asked the question on gis.stackexhange.com. (http://gis.stackexchange.com/questions/191565/vertical-datum-conversion-of-a-raster-with-gdal)

Thanks a lot for your help!!