[Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters
Noel Zinn (cc)
ndzinn at comcast.net
Fri Oct 15 13:57:24 EST 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
>
More information about the Proj
mailing list