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

[Mikael Rittri]

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,

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.

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)