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

[Mikael Rittri]

Noel and Oscar,

I have followed your conversation. Thanks for your interest. 

Noel wrote:

> I've coded this example in Matlab and agree with your final results
> > 66d4'54.705"W   9d34'49.001"N   180.499
> 
> > Nevertheless, I disagree with your first intermediate computation
> >  2464278.090  -5783490.207  974643.507
> 
> My intermediate results are
> =>  2464278.08979296  -5783490.39620226   974642.385933922
>
> These are decimetric differences.  Too much.  Concerned, I computed 
> the same in a commercial offering, Blue Marble Geographic Calculator.  
> Blue Marble agrees with Matlab to the 4 decimal places I happened to 
> have the precision set.  So, I believe that Proj.4 is wrong in this case.  
> You might consider writing a bug report.  (BTW, these differences have 
> nothing to do the use of a single, small-angle rotation matrix versus 
> 3 concatenated single-angle matrices as done by ESRI.) 

I noticed that Oscar got the same results as you. So I tried with 
Carmenta Engine, and I, too, got the same results as you and Oscar, 
apart from a 2.5 mm difference (probably caused by our using three
concatenated single-angle matrices).  

As I started to write a bug report to Proj.4, I thought I would 
run the following code in C++: 

   double dx = 0.0;
   double dy = 0.0;
   double dz = 0.0;
   double arcsecToRadians = 3.141592653589793 / (180. * 60. * 60.);
   double rx = 5.226 * arcsecToRadians;
   double ry = 1.238 * arcsecToRadians;
   double rz = -2.381 * arcsecToRadians;
   double scaleFactor = 1.0 - 5.109e-6;

   double xin =  2464351.59;
   double yin = -5783466.61;
   double zin =   974809.81;

   double xout = scaleFactor * ( 1.0 * xin  - rz * yin  +  ry * zin) + dx;
   double yout = scaleFactor * (  rz * xin  + 1.0 * yin -  rx * zin) + dy;
   double zout = scaleFactor * ( -ry * xin  + rx * yin  + 1.0 * zin) + dz;

which is my (personal) interpretation of the Position Vector Transformation
(as described in EPSG Guidance Note 7.2, section 2.4.3.2.1, pages 109-110.)

To my surprise, now I got 

	xout:  2464278.0897929627	
	yout: -5783490.2071627714	
	zout:   974643.50748968683	

which agrees with cs2cs, but not with Noel/Oscar/Carmenta Engine ...
So I had to cancel my bug report.

Now I am just completely confused. Does my C++ code above contain
some misunderstanding of the Position Vector Transformation,
a misunderstanding that is also present in the cs2cs implementation? 
Or is cs2cs correct after all? 

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

-----Original Message-----
From: proj-bounces at lists.maptools.org [mailto:proj-bounces at lists.maptools.org] On Behalf Of Noel Zinn (cc)
Sent: den 15 oktober 2010 20:57
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters

[ see http://lists.maptools.org/pipermail/proj/2010-October/005441.html ]