[Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters
Mikael.Rittri at carmenta.com
Tue Oct 12 03:45:48 EST 2010
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
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
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
To test it, we use the test point given in section 18.104.22.168,
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.
More information about the Proj