[Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters
Mikael.Rittri at carmenta.com
Mon Oct 18 07:50:30 EST 2010
Oops. Noel wrote:
> Your C++ code is correct.
> The value of rX is inconsistent in your documentation, changing from 5.266 to 5.226 arcseconds.
> This would explain why it affects Y and Z but not X.
Yes indeed. In the list of parameters in the beginning of my example,
I gave the correct value of (minus) 5.266, but when I used this rotation in cs2cs commands,
I used the wrong value of 5.226 instead. Noel and Oscar have used the correct value,
which I also happened to do (without noticing) when I tried the computation in Carmenta
Engine. But in my C++ code, I reverted to the wrong value again.
I'd like to correct my example:
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.266,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.396 974642.386
ii) Now subtract the output row from the input row, and we get
73.500 23.786 167.424
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.385 -192.802
The last three values are the question marks in our 7-parameter
To test it, we use the test point given in section 188.8.131.52,
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.385,-192.802,5.266,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.514
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.
In my first, wrongly computed example, I had a mismatch of about 1.1 cm.
I think it is gone now. Since I used the wrong rx rotation twice, I think
the two errors nearly, but not quite, cancelled each other.
Of course, given the number of decimals used in GN 7.2 for the correct
result, there might still be a mismatch up to 0.001 arc seconds in long/lat
and up to 1 cm in height.
Noel also wrote:
> Is it something going on under the hood in Proj.4,
> like an embedded call to geodetic <=> geocentric for some reason?
I am not sure, but I think so. But I think this would cause less
than 1 millimeter error, at least in my example.
From: Noel Zinn (cc) [mailto:ndzinn at comcast.net]
Sent: den 18 oktober 2010 13:54
To: Mikael Rittri
Subject: Offline - Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters
Your C++ code is correct.
The value of rX is inconsistent in your documentation, changing from 5.266 to 5.226 arcseconds.
This would explain why it affects Y and Z but not X.
Noel Zinn, Principal, Hydrometronics LLC
+1-832-539-1472 (office), +1-281-221-0051 (cell)
noel.zinn at hydrometronics.com (email)
More information about the Proj