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

Mikael Rittri Mikael.Rittri at carmenta.com
Mon Oct 18 06:06:03 EST 2010

```Noel and Oscar,

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 ]
```