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

Fri Oct 15 11:57:24 PDT 2010

Mikael,

Your provision of a worked example "compressing" a M-B transformation into a 
7-parameter transformation is a real contribution.  Thanks ... and thanks to 
Melita Kennedy, too.

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.)

We agree in the final result because the next intermediate steps (computing 
dX, dY, dZ, which are also different in my case) are fortuitously self 
correcting (at least at the single evaluation point).

So, thanks again, Mikael, for this useful procedure.  More comments later.

Regards,
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)

--------------------------------------------------
From: "Mikael Rittri" <Mikael.Rittri at carmenta.com>
Sent: Tuesday, October 12, 2010 3:45 AM
To: "PROJ.4 and general Projections Discussions" <proj at lists.maptools.org>
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
> 