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

[OvV_HN]

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
Sun Oct 17 07:44:35 EST 2010

.....

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).

.......


> 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

---------------------------------

REPLY:

We do agree. Something is somewhere wrong. But what and where?
If I were to get a bug report stating that in a procedure spanning several 
operations the result departs from the expected value, I would throw this 
report away.
Every single step should be tested against independent results from others, 
not just the end result of a concatenation of a handful of steps.

Several steps are clear, for instance:
- The right choice of rotation model is made.
- BTW: rotation order has no effect if you work with approximations in the 
rotation matrices, like PROJ does.
- The datum transformation functions of PROJ seem to be correct.
- At least the core of the geocentric->geodetic conversion function is OK.

But where does it go wrong?
- Are all parameters passed correctly? Is somewhere a lot of precision lost 
because at a higher level not enough decimals are set?
- Are the global structs filled correctly by initialization functions?
- Is the software used as it should be? Take cs2cs: you can do incredible 
things with it, but next to nothing is documented (see its man page).

Oscar van Vlijmen