[Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters
Mikael Rittri
Mikael.Rittri at carmenta.com
Mon Oct 18 06:06:03 EST 2010
Noel and Oscar,
I have followed your conversation. Thanks for your interest.
Noel wrote:
> 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.)
I noticed that Oscar got the same results as you. So I tried with
Carmenta Engine, and I, too, got the same results as you and Oscar,
apart from a 2.5 mm difference (probably caused by our using three
concatenated single-angle matrices).
As I started to write a bug report to Proj.4, I thought I would
run the following code in C++:
double dx = 0.0;
double dy = 0.0;
double dz = 0.0;
double arcsecToRadians = 3.141592653589793 / (180. * 60. * 60.);
double rx = 5.226 * arcsecToRadians;
double ry = 1.238 * arcsecToRadians;
double rz = -2.381 * arcsecToRadians;
double scaleFactor = 1.0 - 5.109e-6;
double xin = 2464351.59;
double yin = -5783466.61;
double zin = 974809.81;
double xout = scaleFactor * ( 1.0 * xin - rz * yin + ry * zin) + dx;
double yout = scaleFactor * ( rz * xin + 1.0 * yin - rx * zin) + dy;
double zout = scaleFactor * ( -ry * xin + rx * yin + 1.0 * zin) + dz;
which is my (personal) interpretation of the Position Vector Transformation
(as described in EPSG Guidance Note 7.2, section 2.4.3.2.1, pages 109-110.)
To my surprise, now I got
xout: 2464278.0897929627
yout: -5783490.2071627714
zout: 974643.50748968683
which agrees with cs2cs, but not with Noel/Oscar/Carmenta Engine ...
So I had to cancel my bug report.
Now I am just completely confused. Does my C++ code above contain
some misunderstanding of the Position Vector Transformation,
a misunderstanding that is also present in the cs2cs implementation?
Or is cs2cs correct after all?
Mikael Rittri
Carmenta AB
Sweden
http://www.carmenta.com
-----Original Message-----
From: proj-bounces at lists.maptools.org [mailto:proj-bounces at lists.maptools.org] On Behalf Of Noel Zinn (cc)
Sent: den 15 oktober 2010 20:57
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters
[ see http://lists.maptools.org/pipermail/proj/2010-October/005441.html ]
More information about the Proj
mailing list