[Proj] Geoid models egm96 et egm08 in PROJ4

[Alain Orsoni]

Hi,

We think there could be a problem in the egm96_15.gtx and egm08_25.gtx files provided on ftp://ftp.remotesensing.org/proj/vdatum/

We noticed a difference between the values we used to get and the values provided by PROJ4 using these gtx files.
It seems that both files have a shift (although it is different for each of the files.

Probably the reason is related to the egm08togtx.py and egm96togtx.py :

=============   in egm96togtx.py the line
ds_gtx.SetGeoTransform( (-180.25,0.25,0,90.125,0,-0.25) )
should probably become
ds_gtx.SetGeoTransform( (-180.125,0.25,0,90.125,0,-0.25) )
assuming that the grid value is given for the center of a cell (0.25°x0.25° cell). 

This hypothesis seems to be confirmed by the following verification:
on the web site http://www.unavco.org/community_science/science-support/geoid/
which is a geoid height calculator, if you enter 0 0 and 0 values for long lat and height you will get a -17.162 orthometric height

If you use the following command

cs2cs +init=epsg:4326 +to +proj=longlat +datum=WGS84 +geoidgrids=egm96_15.gtx +units=m +no_defs
you will get the following result
0 0 0
0dE	0dN -17.120
while if you shift your value of -0.125
-0.125 0 0
0d7'30"W	0dN -17.160


=============  in egm08togtx.py the line
    ps_25 = 2.5 / 60.0
    ds_gtx.SetGeoTransform( (-180 - ps_25,ps_25,0,90 + ps_25,0,-1 * ps_25) )

should become

    ds_gtx.SetGeoTransform( (-180 - (ps_25/2),ps_25,0,90 + (ps_25/2),0,-1 * ps_25) )

=================

If some one as encountered the same problem and confirm or could unconfirm this hypothesis, it would be helpful.

Alain ORSONI