Huge differences between results for 7-param datum transform

Question

I'm trying to transform a WGS84 lat/lon point

11d42'32.10629" E 5d12'56.75371" S

in three different software packages (Proj4, GeoTrans and Leica GeoOffice) but I receive very big differences between the results - around few hundred meters in X and Y! Note that these variations appears only with 7 parameters not with 3.

Proj4

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=utm +zone=32 +ellps=clrk80 +towgs84=178.3,316.7,131.5,-5.278,-6.077,-10.9792,-19.166 +south +units=m +no_defs

---

GeoTrans

Delta X 178.3 
Delta Y 316.7 
Delta Z 131.5 
Rot X -5.278 
Rot Y -6.077 
Rot Z -10.9792 
SF = -19.166 / 0.999980834 (0.000019166)

---

Leica GeoOffice

---

The results are respectively:

1. 800392 9422525
2. 800306 9422840
3. 800941 9422891

I believe all the 3 packages use same mathematical methods for 7-param transform (Bursa-Wolf method). What could be the cause of this huge variation?

Answer 1

First, Proj4 uses what EPSG calls the "Position Vector" version of the 7 parameter method. It's possible that GeoTrans and Leica GeoOffice use the other version which EPSG called "Coordinate Frame". Both methods are equivalent, but the rotation matrices are different and the signs of the angular parameters have to be changed.

Second, thank you for providing the screen shot of the transformation definition in Leica GeoOffice. Either the Proj4 definition or this definition is defined incorrectly. The defining parameters Ellipsoid A and Ellipsoid B should be switched. Currently, this transformation converts from WGS84 to Congo60. In Proj4, note that the option is +towgs84, so it's defined as FROM Congo60 TO WGS84. To change the direction in the definition, change the signs of ALL the parameters. Also check the help for GeoOffice and see whether the SF parameter wants the parts per million version or the already converted-to-a-scale-factor value.

I don't know about GeoTrans--do you mean the NGA software? Anyway, hopefully you can get GeoOffice and Proj4 to match.

Answer 2

For proj4, the signs of the parameters have to be reversed.

See this definition page:

http://www.spatialreference.org/ref/epsg/62826405/prettywkt/

GEOGCS["Pointe Noire (deg)",
    DATUM["Congo 1960 Pointe Noire",
        SPHEROID["Clarke 1880 (IGN)",6378249.2,293.4660212936269,
            AUTHORITY["EPSG","7011"]],
        TOWGS84[-178.3,-316.7,-131.5,5.278,6.077,10.979,3.953271276531849],
        AUTHORITY["EPSG","6282"]],
    PRIMEM["Greenwich",0.0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.017453292519943295],
    AXIS["Geodetic latitude",NORTH],
    AXIS["Geodetic longitude",EAST],
    AUTHORITY["EPSG","62826405"]]

I'm not sure where you got your values from.

QGIS defines EPSG:28232 Point Noire UTM 32S as:

+proj=utm +zone=32 +south +a=6378249.2 +b=6356515 +towgs84=-148,51,-291,0,0,0,0 +units=m +no_defs

and the 7 parameters result should not be far away from that.

---

EDIT

With GDAL 1.10, I get the following results:

4326-proj-3 parameters:
cs2cs +init=epsg:4326 +to +proj=utm +zone=32 +south +a=6378249.2 +b=6356515 +towgs84=-148,51,-291,0,0,0,0 +units=m +no_defs
800232.21   9423131.96 -1.76
4326-proj-7 parameters from spatialrefrence.org:
cs2cs +init=epsg:4326 +to +proj=utm +zone=32 +south +a=6378249.2 +b=6356515 +towgs84=-178.3,-316.7,-131.5,5.278,6.077,10.979,3.953271276531849 +units=m +no_defs
800230.13   9423133.46 91.31
4326-proj-7 parameters from proj4 datum_shift.csv (EPSG:1802):
cs2cs +init=epsg:4326 +to +proj=utm +zone=32 +south +a=6378249.2 +b=6356515 +towgs84=-178.3,-316.7,-131.5,5.278,6.077,10.979,19.166 +units=m +no_defs
800230.13   9423133.40 -5.72

So less than 3 metres apart. Note that the E and S must follow the degrees without blanks.

And, just for the record: setting towgs84 to zero results in:

800310.94   9422829.37 -109.32

while transforming from Point Noire degrees EPSG:4262 to Point Noire UTM results in:

800311.21   9422892.49 0.00