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

[Mikael Rittri]

A Molodensky-Badekas datum shift is represented by 10 parameters. 
But it can be compressed to a 7-parameter representation. 

As Noel Zinn has explained (http://lists.maptools.org/pipermail/proj/2010-October/005424.html),
the 10-parameter representation has advantages over the 7-parameter 
one, for people who derive, analyze and compare datum shifts. 

But for people who just use a datum shift to transform coordinates,
the 7-parameter representation gives the same results. In this 
limited sense, the compression is lossless, which is useful if 
your software (like Proj.4) does not support 10-parameter M-B.

I learned how to do the compression from Melita Kennedy (thanks 
again!).  And I have now realized that one can persuade Proj.4 to 
do most of the work of the compression process.  Since I thought 
this is rather nice, I'd like to explain it by an example. 

Example: M-B datum shift EPSG:1096, "La Canoa to WGS 84 (2)",
which has the same parameters as EPSG:1771, "La Canoa to REGVEN (1)".
    
    Source ellipsoid: International 1924.
    Target ellipsoid: WGS 84.
	
    dX = -270.933 m
    dY =  115.599 m
    dZ = -360.226 m
    RX = -5.266 arc sec 
    RY = -1.238 arc sec 
    RZ =  2.381 arc sec 
    dS = -5.109 ppm
    eX =  2464351.59 m
    eY = -5783466.61 m
    eZ =   974809.81 m

where eX,eY,eZ are the geocentric coordinates of the so-called 
evaluation point.

In principle, the equivalent 7-parameter datum shift should have 
the same rotations RX, RY, RZ and the same scale difference dS.
But there is a catch: the sign convention for the rotations. 
M-B uses the same sign convention as Coordinate Frame Rotation, 
while cs2cs uses the opposite convention of the Position Vector 
Transform, so we must negate the three rotations.  So, in cs2cs 
syntax, the 7-parameter datum shift we seek has the form

    +towgs84=?,?,?,5.266,1.238,-2.381,-5.109
	
To find the question marks, we have to do three steps.  

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.226,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.207  974643.507

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

        73.500        23.597     166.303

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.196    -193.923
   
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.196,-193.923,5.226,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.499

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 - apart from the 
unexplained difference of about 1.1 cm in the ellipoidal height.  
  
Regards,
Mikael Rittri
Carmenta AB
Sweden
www.carmenta.com