[Proj] Geoid models egm96 et egm08 in PROJ4

[Charles Karney]

Corrected versions of egm96_15.gtx and egm08_25.gtx are available at

   https://sarnoff.sharefile.com/d/sb23d64908da48ffb

These were produced directly from the NGA spherical harmonic expansions
using the GeoidToGTX program (a sample program distributed with
GeographicLib).  These files correct the following issues with the PROJ4
files:

(1) +/- 5mm quantization errors in egm96.
(2) incorrect origins in egm96 and egm08

GeoidToGTX faithfully reproduces the grids given by the harmonic
synthesis programs (in Fortran) provided by NGA.

(The link above should be valid "indefinitely".  However, I can't be
sure that it will work beyond about a year.)

   --Charles

On 2014-02-12 08:29, Alain Orsoni wrote:
> 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