12

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

Screenshot


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?

  • In Geo Office, is the Clarke 1880 IGN ellipsoid a custom definition? Can you post its parameters whether it is or isn't? – mkennedy Sep 26 '13 at 19:29
  • I've received the Leica screenshot from other people and don't really know yet. – WindRider Sep 28 '13 at 13:59
9

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.

  • 1
    The Y (northing/latitude) value could be because the Clarke 1880 values aren't quite the same, but 60 m seems too large. – mkennedy Sep 23 '13 at 19:45
  • 1
    @mkennedy: can you do the same calculations I did in Arcgis to see if we can align without Leica? – AndreJ Sep 24 '13 at 19:10
  • 1
    @AndreJoost, sure. Using the Esri versions of EPSG 1802 and 28232, I got easting: 800230.139 northing:9423133.413. – mkennedy Sep 24 '13 at 19:22
  • 1
    I tried a several variants to see if I could reproduce the Leica result and wasn't able to. Same issue that you found: Y off by 60+ m. – mkennedy Sep 26 '13 at 19:15
  • 1
    I am confident that Proj4 (and ArcGIS!) solutions are correct. The Proj4 source is open. Although the ArcGIS code isn't open, I've dealt with it for about 15 years. I also added a custom transformation to Geotrans 3.0 and got the same result as Proj4/ArcGIS. – mkennedy Sep 26 '13 at 19:40
4

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
  • As you can see, with same parameters there is a huge difference (>600 m!) from Leica. I need to know why this happens and how to obtain the same results. Could it be from the method of conversion from/to geocentric coordinates used in Proj4? It has very complicated formulas. Maybe Leica use different method? – WindRider Sep 23 '13 at 16:05
  • 1
    Do Leica offer 3-parameter-conversion? And what do they compute for it? – AndreJ Sep 23 '13 at 17:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.