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

Tue Oct 12 07:37:26 PDT 2010

The "evaluation point" is the Datum Origin, a rather important item.  In regard to the "unexplained" difference in ellipsoid height results, my guess is the approximations made by the inverse algorithm from Geocentric coordinates.

Personally, I prefer the elegance of incorporating the Datum Origin.  

I have also been told by Dr. Muneendra Kumar that Badekas had nothing to do with Molodensky's original method.  I believe they were in school together at the same time, and Prof. Molodensky preceeded them on another continent by some years if not decades.

Clifford J. Mugnier, C.P., C.M.S.
Chief of Geodesy,
Center for GeoInformatics
Department of Civil Engineering 
Patrick F. Taylor Hall 3223A
LOUISIANA STATE UNIVERSITY 
Baton Rouge, LA  70803
Voice and Facsimile:  (225) 578-8536 [Academic] 
Voice and Facsimile:  (225) 578-4578 [Research] 
Cell: (225) 238-8975 [Academic & Research]
Honorary Life Member of the 
Louisiana Society of Professional Surveyors 
Fellow Emeritus of the ASPRS 
Member of the Americas Petroleum Survey Group

________________________________

From: proj-bounces at lists.maptools.org on behalf of Mikael Rittri
Sent: Tue 12-Oct-10 03:45
To: PROJ.4 and general Projections Discussions
Subject: [Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters

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
_______________________________________________
Proj mailing list
Proj at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20101012/9d6198fe/attachment.html>