[Proj] Geoid models egm96 et egm08 in PROJ4

[Mikael Rittri]

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

I agree with you, and I have reported this problem earlier, as 
 
   https://trac.osgeo.org/proj/ticket/126

Charles Karney has also confirmed it, in 

   http://lists.maptools.org/pipermail/proj/2012-April/006273.html


Best regards,

Mikael Rittri
Carmenta
Sweden
http://www.carmenta.com

-----Original Message-----
From: proj-bounces at lists.maptools.org [mailto:proj-bounces at lists.maptools.org] On Behalf Of Alain Orsoni
Sent: Wednesday, February 12, 2014 2:30 PM
To: proj at lists.maptools.org
Subject: [Proj] Geoid models egm96 et egm08 in PROJ4

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

_______________________________________________
Proj mailing list
Proj at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj