[Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekasdatum shift into 7 parameters

[Noel Zinn (cc)]

Thanks, Oscar, for the useful tour of geocentric <=> geodetic conversions in 
Proj.4.  The geodetic => geocentric piece is (should be) exact in any 
formulation.  The geocentric => geodetic piece lends itself to an iterative 
solution for precision if you find yourself in a place humans can't inhabit 
normally.  Some of this ground we've covered previously in this forum, 
notably:

http://lists.maptools.org/pipermail/proj/2009-January/004315.html

To help me understand whether we agree of not, let me decompose Mikael 
Rittri's problem simply.

A 7-parameter datum transformation (geodetic <=> geodetic) has 3 parts: (1) 
geodetic => geocentric, (2) matrix manipulation on 3D geocentric (ECEF) 
coordinates involving translation, rotation and scale change, and then (3) 
geocentric => geodetic (which is the only inexact piece of the process, but 
you can get any precision you desire).

Now, Mikael's first step of his (Melita's) compression process is just (2) 
of the above, the geocentric => geocentric matrix manipulation in 3D 
Cartesian space using the following Proj.4 command (if I understand 
correctly):

>cs2cs -f "%.3f" +proj=geocent +towgs84=0,0,0,5.226,1.238,-2.381,-5.109 +to 
>+datum=WGS84 +proj=geocent

His input is  2464351.59   -5783466.61   974809.81

His output is  2464278.090  -5783490.207  974643.507

To reproduce this in Matlab I pulled out parts (1) and (3) from my 
7-parameter transformation code and using Mikael's input I get a different 
output (2464278.090  -5783490.396   974642.386) which is the same as that 
served up by your own Molodensky-Badekas code:

> Let's see. In the problem at hand my Molodenksy-Badekas routine arrives at
> an X,Y,Z of:
> X = 2464278.0897929626; Y = -5783490.396202258; Z = 974642.3859339222;

So, we agree then?  The Proj.4 geocentric => geocentric routine (cs2cs -f 
"%.3f" +proj=geocent +towgs84=0,0,0,5.226,1.238,-2.381,-5.109 +to 
+datum=WGS84 +proj=geocent) has a decimetric difference with your code, my 
code and with Blue Marble.

I have verified that this difference is not due to the difference between 
the single rotation matrix normally used (small angle approximation) and the 
3 concatenated rotation matrices used by ESRI (for example), which is only 
millimetric in quantity.  Is it something going on under the hood in Proj.4, 
like an embedded call to geodetic <=> geocentric for some reason?  I don't 
know.  We could debate whether 19 centimeters in Y and 12 centimeters in Z 
is significant or not.  For a straight-up matrix manipulation, I think it 
is.

Regards,
Noel

Noel Zinn, Principal, Hydrometronics LLC
+1-832-539-1472 (office), +1-281-221-0051 (cell)
noel.zinn at hydrometronics.com (email)
http://www.hydrometronics.com (website)

--------------------------------------------------
From: "OvV_HN" <ovv at hetnet.nl>
Sent: Sunday, October 17, 2010 3:15 AM
To: "PROJ.4 and general Projections Discussions" <proj at lists.maptools.org>
Subject: [Proj] Using Proj.4 to compress a 10-parameter 
Molodensky-Badekasdatum shift into 7 parameters

> IN REPLY TO:
>
> [Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekasdatum 
> shift
> into 7 parameters
> Noel Zinn (cc) ndzinn at comcast.net
> Sat Oct 16 09:42:42 EST 2010
>
> .......
>
>> I haven't tested the geocentric -> geodetic code from PROJ4's geocent.c,
>> but
>> it certainly looks good.
>
> .......
>
> I contend that the correct output (delivered by my Matlab scripts and by
> Blue Marble) is
> =>  2464278.090  -5783490.396   974642.386 (to the same precision)
>
> If my assertion is correct, then Proj.4 has a slight problem.  It's not
> blame, just a request for an independent verification of this particular
> step (i) and (if verified) another possible improvement for Proj.4.
>
> Noel
>
>
>
> REPLY:
>
> I do get exactly the same output as stated above, so up to that point 
> there
> is nothing wrong I guess.
> Let's look at PROJ4's geocent. I took the code apart and tested both
> geocentric->geodetic functions.
> Indeed, PROJ4's geocent.c has two functions. Which one is actually used,
> depends on a compiler directive.
> The first code is from Ralph Toms. This code is optimized for speed. It 
> uses
> a small fudge factor (AD_C). This introduces a small error, but it 
> enhances
> the error stability especially for extremely high altitudes.
> The second code is an iterative procedure, which gives nearly exact 
> results
> (the iteration cut-off error is 1e-12).
>
> Let's see. In the problem at hand my Molodenksy-Badekas routine arrives at
> an X,Y,Z of:
> X = 2464278.0897929626; Y = -5783490.396202258; Z = 974642.3859339222;
> Converted to lat, lon, height on WGS84 with the exact Borkowski procedure
> gives:
> 9.580277979760344, -66.08186260475033, 180.5140539117856
> Converted with the iterative procedure from geocent.c gives in my test:
> 9.580277979760346, -66.08186260475033, 180.514053911902
> Converted with the Toms procedure from geocent.c gives in my test:
> 9.580277981233836, -66.08186260475033, 180.51408141944557
> The latter test differs by 5.3 micro-arcseconds in the latitude with the
> Borkowski method and 27.5 micrometers in the height.
>
> I believe PROJ4's geocent uses by default the somewhat inaccurate Toms
> method for converting geocentric to geodetic.
> But at least in my test it seems that the error is rather insignificantly
> small.
> So the core of geocent.c seems OK.
> Are the global structs filled correctly and are all parameters passed
> correctly?
>
>
> Oscar van Vlijmen
>
>
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
>