[Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekasdatum shift into 7 parameters
Noel Zinn (cc)
ndzinn at comcast.net
Sun Oct 17 07:44:35 EST 2010
Thanks, Oscar, for the useful tour of geocentric <=> geodetic conversions in
Proj.4. The geodetic => geocentric piece is (should be) exact in any
formulation. The geocentric => geodetic piece lends itself to an iterative
solution for precision if you find yourself in a place humans can't inhabit
normally. Some of this ground we've covered previously in this forum,
notably:
http://lists.maptools.org/pipermail/proj/2009-January/004315.html
To help me understand whether we agree of not, let me decompose Mikael
Rittri's problem simply.
A 7-parameter datum transformation (geodetic <=> geodetic) has 3 parts: (1)
geodetic => geocentric, (2) matrix manipulation on 3D geocentric (ECEF)
coordinates involving translation, rotation and scale change, and then (3)
geocentric => geodetic (which is the only inexact piece of the process, but
you can get any precision you desire).
Now, Mikael's first step of his (Melita's) compression process is just (2)
of the above, the geocentric => geocentric matrix manipulation in 3D
Cartesian space using the following Proj.4 command (if I understand
correctly):
>cs2cs -f "%.3f" +proj=geocent +towgs84=0,0,0,5.226,1.238,-2.381,-5.109 +to
>+datum=WGS84 +proj=geocent
His input is 2464351.59 -5783466.61 974809.81
His output is 2464278.090 -5783490.207 974643.507
To reproduce this in Matlab I pulled out parts (1) and (3) from my
7-parameter transformation code and using Mikael's input I get a different
output (2464278.090 -5783490.396 974642.386) which is the same as that
served up by your own Molodensky-Badekas code:
> Let's see. In the problem at hand my Molodenksy-Badekas routine arrives at
> an X,Y,Z of:
> X = 2464278.0897929626; Y = -5783490.396202258; Z = 974642.3859339222;
So, we agree then? The Proj.4 geocentric => geocentric routine (cs2cs -f
"%.3f" +proj=geocent +towgs84=0,0,0,5.226,1.238,-2.381,-5.109 +to
+datum=WGS84 +proj=geocent) has a decimetric difference with your code, my
code and with Blue Marble.
I have verified that this difference is not due to the difference between
the single rotation matrix normally used (small angle approximation) and the
3 concatenated rotation matrices used by ESRI (for example), which is only
millimetric in quantity. Is it something going on under the hood in Proj.4,
like an embedded call to geodetic <=> geocentric for some reason? I don't
know. We could debate whether 19 centimeters in Y and 12 centimeters in Z
is significant or not. For a straight-up matrix manipulation, I think it
is.
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: "OvV_HN" <ovv at hetnet.nl>
Sent: Sunday, October 17, 2010 3:15 AM
To: "PROJ.4 and general Projections Discussions" <proj at lists.maptools.org>
Subject: [Proj] Using Proj.4 to compress a 10-parameter
Molodensky-Badekasdatum shift into 7 parameters
> IN REPLY TO:
>
> [Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekasdatum
> shift
> into 7 parameters
> Noel Zinn (cc) ndzinn at comcast.net
> Sat Oct 16 09:42:42 EST 2010
>
> .......
>
>> I haven't tested the geocentric -> geodetic code from PROJ4's geocent.c,
>> but
>> it certainly looks good.
>
> .......
>
> I contend that the correct output (delivered by my Matlab scripts and by
> Blue Marble) is
> => 2464278.090 -5783490.396 974642.386 (to the same precision)
>
> If my assertion is correct, then Proj.4 has a slight problem. It's not
> blame, just a request for an independent verification of this particular
> step (i) and (if verified) another possible improvement for Proj.4.
>
> Noel
>
>
>
> REPLY:
>
> I do get exactly the same output as stated above, so up to that point
> there
> is nothing wrong I guess.
> Let's look at PROJ4's geocent. I took the code apart and tested both
> geocentric->geodetic functions.
> Indeed, PROJ4's geocent.c has two functions. Which one is actually used,
> depends on a compiler directive.
> The first code is from Ralph Toms. This code is optimized for speed. It
> uses
> a small fudge factor (AD_C). This introduces a small error, but it
> enhances
> the error stability especially for extremely high altitudes.
> The second code is an iterative procedure, which gives nearly exact
> results
> (the iteration cut-off error is 1e-12).
>
> Let's see. In the problem at hand my Molodenksy-Badekas routine arrives at
> an X,Y,Z of:
> X = 2464278.0897929626; Y = -5783490.396202258; Z = 974642.3859339222;
> Converted to lat, lon, height on WGS84 with the exact Borkowski procedure
> gives:
> 9.580277979760344, -66.08186260475033, 180.5140539117856
> Converted with the iterative procedure from geocent.c gives in my test:
> 9.580277979760346, -66.08186260475033, 180.514053911902
> Converted with the Toms procedure from geocent.c gives in my test:
> 9.580277981233836, -66.08186260475033, 180.51408141944557
> The latter test differs by 5.3 micro-arcseconds in the latitude with the
> Borkowski method and 27.5 micrometers in the height.
>
> I believe PROJ4's geocent uses by default the somewhat inaccurate Toms
> method for converting geocentric to geodetic.
> But at least in my test it seems that the error is rather insignificantly
> small.
> So the core of geocent.c seems OK.
> Are the global structs filled correctly and are all parameters passed
> correctly?
>
>
> Oscar van Vlijmen
>
>
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
>
More information about the Proj
mailing list