# [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 07:50:30 EST 2010

```Oops. Noel wrote:

> Your C++ code is correct.
> The value of rX is inconsistent in your documentation, changing from 5.266 to 5.226 arcseconds.
> This would explain why it affects Y and Z but not X.

Yes indeed.  In the list of parameters in the beginning of my example,

http://lists.maptools.org/pipermail/proj/2010-October/005436.html

I gave the correct value of (minus) 5.266, but when I used this rotation in cs2cs commands,
I used the wrong value of 5.226 instead.  Noel and Oscar have used the correct value,
which I also happened to do (without noticing) when I tried the computation in Carmenta
Engine.  But in my C++ code, I reverted to the wrong value again.

Sorry.

I'd like to correct my example:

<CORRECTED EXAMPLE>

...

i) Make a 7-parameter "something" by temporarily replacing
the question marks by zeros.  It is not really a datum
shift, but we might call it a 7-parameter compressor.
Use it in a cs2cs command like this

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

As input, use the evaluation point eX,eY,eZ of the M-B datum shift:

2464351.59   -5783466.61   974809.81

We get the output

2464278.090  -5783490.396  974642.386

ii) Now subtract the output row from the input row, and we get

73.500        23.786     167.424

iii) Then add the original dX,dY,dZ from the M-B datum shift,
which were -270.933,115.599,-360.226, and we get

-197.433       139.385    -192.802

The last three values are the question marks in our 7-parameter
datum shift.

To test it, we use the test point given in section 2.4.4.1,
pages 114-115, EPSG Guidance Note 7.2, http://www.epsg.org/guides/docs/G7-2.pdf

>cs2cs +ellps=intl +proj=longlat +towgs84=-197.433,139.385,-192.802,5.266,1.238,-2.381,-5.109 +to +datum=WGS84 +proj=longlat
66d4'48.091"W   9d35'0.386"N    201.46
66d4'54.705"W   9d34'49.001"N   180.514

According to GN 7.2, the result should have been

66°04'54.705"W	9°34'49.001"N   180.51 m

so we have achieved a lossless compression.
</CORRECTED EXAMPLE>

In my first, wrongly computed example, I had a mismatch of about 1.1 cm.
I think it is gone now. Since I used the wrong rx rotation twice, I think
the two errors nearly, but not quite, cancelled each other.
Of course, given the number of decimals used in GN 7.2 for the correct
result, there might still be a mismatch up to 0.001 arc seconds in long/lat
and up to 1 cm in height.

Noel also wrote:
> Is it something going on under the hood in Proj.4,
> like an embedded call to geodetic <=> geocentric for some reason?

I am not sure, but I think so. But I think this would cause less
than 1 millimeter error, at least in my example.

Best regards,
Mikael Rittri
Carmenta AB
http://www.carmenta.com

-----Original Message-----
From: Noel Zinn (cc) [mailto:ndzinn at comcast.net]
Sent: den 18 oktober 2010 13:54
To: Mikael Rittri
Subject: Offline - Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters

Mikael,

The value of rX is inconsistent in your documentation, changing from 5.266 to 5.226 arcseconds.

This would explain why it affects Y and Z but not X.

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)

```