Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Need conversion from Tanzanian ARC 1960 to WGS84 and vice-versa #1110

Closed
gnusupport opened this issue Sep 5, 2018 · 23 comments
Closed

Need conversion from Tanzanian ARC 1960 to WGS84 and vice-versa #1110

gnusupport opened this issue Sep 5, 2018 · 23 comments

Comments

@gnusupport
Copy link

Hello,

I need conversion from Tanzanian ARC1960 geodetic system to WGS84, and vice versa.

Please contact us over http://www.startyourowngoldmine.com and arrange your donation, it should be something reasonable.

Thank you.

@gnusupport
Copy link
Author

And conversion shall be truthful, I think that GEOTRANS software is doing it well, but that is not good API option.

@kbevers
Copy link
Member

kbevers commented Sep 5, 2018

You need to help us help you. Please supply a system definition of ARC1960 - otherwise this will get nowhere.

Also, you probably have a much better chance of getting help if you don't ask for money via a dodgy looking website.

@gnusupport
Copy link
Author

@kbevers -- you contact me, and propose donation from my side to you. There is misunderstanding on your side, so let us be professional.

Can this link help you: http://georepository.com/datum_6210/Arc-1960.html ?

@kbevers
Copy link
Member

kbevers commented Sep 5, 2018

I'm sorry I misunderstood you. Your question has been asked on the mailing list previously. See this thread: http://lists.maptools.org/pipermail/proj/2009-August/004849.html

@gnusupport
Copy link
Author

Yes, thank you. Do you think that the link you have given may be fully applied by proj programs?

Or better: can I use proj to make that transformation without external programs?

@kbevers
Copy link
Member

kbevers commented Sep 6, 2018

Yes, thank you. Do you think that the link you have given may be fully applied by proj programs?

Yes.

Or better: can I use proj to make that transformation without external programs?

Also yes.

The transformation doesn't look particularly difficult but without test coordinates it is not possible to verify if the transformation is configured correctly. Can you provide coordinates for a few points in both ARC 1960 and WGS84?

@gnusupport
Copy link
Author

Thank you much. Here it is ARC 1960 for Tanzania (which is not quite same as ARC 1960 for Kenya or Uganda):

Latitude, longitude:
6° 55' 57" S and 37° 23' 30" E
6° 55' 57" S and 37° 23' 36" E

and those both are translated by GEOTRANS to:

6° 56' 4.7" S and 37° 23' 32.9" E
6° 56' 4.7" S and 37° 23' 38.9" E

and I can just believe that GEOTRANS is doing good job.

@kbevers
Copy link
Member

kbevers commented Sep 7, 2018

I have made a short gie file to test the transformation based on your supplied coordinates. It's displayed below. From the file you can see how to set up the transformation. Beware that the accuracy of the transformation is around 2 m (which is to be expected from an old datum in combination with the molodensky transform). I am closing this ticket since the problem is solved.

<gie>

operation   +proj=molodensky +ellps=clrk80
                             +da=-112.145 +df=-0.54750714e-4
                             +dx=-175 +dy=-23 +dz=-303
tolerance   2 m

accept      37.3916666667   -6.9325000000
expect      37.3924722222   -6.9346388889

accept      37.3933333333   -6.9325000000
expect      37.3941388889   -6.9346388889

</gie>

@kbevers kbevers closed this as completed Sep 7, 2018
@gnusupport
Copy link
Author

cs2cs +proj=molodensky +ellps=clrk80 +da=-112.145 +df=-0.54750714e-4 +dx=-175 +dy=-23 +dz=-303 +to +proj=latlon +datum=WGS84 < file

where file is:

37.3916666667   -6.9325000000                                                       
37.3933333333   -6.9325000000

please help me, what I am doing wrong?

As I am getting:

cs2cs +proj=molodensky +ellps=clrk80 +da=-112.145 +df=-0.54750714e-4 +dx=-175 +dy=-23 +dz=-303 +to +proj=latlon +datum=WGS84 < file
Rel. 5.1.0, June 1st, 2018
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): latitude or longitude exceeded limits
*       * 0.000                                                       
Rel. 5.1.0, June 1st, 2018
<cs2cs>: while processing file: <stdin>, line 2
pj_transform(): latitude or longitude exceeded limits
*       * 0.000                                                       

@kbevers
Copy link
Member

kbevers commented Sep 7, 2018

Use cct instead:

cct +proj=molodensky +ellps=clrk80 +da=-112.145 +df=-0.54750714e-4 +dx=-175 +dy=-23 +dz=-303 file

You will have to add height and time columns to your file or add -z0 -t0 flags to the cct call. Consult the docs for more details.

@gnusupport
Copy link
Author

That is great help, thank you much. It solves so many problems for people.

@gnusupport
Copy link
Author

This location: 6° 55' 57" S, 37° 23' 30" E is -6.932500, 37.391667
and if I put in a file:

-6.932500 37.391667

result from command:

cct -z1000 -t0 +proj=molodensky +ellps=clrk80 +da=-112.145 +df=-0.54750714e-4 +dx=-175 +dy=-23 +dz=-303 +to +proj=latlon +datum=WGS84 file2

is

-6.9329962797   37.3904326033      680.1305

but if I put following notation in the file:

6d55'57"S 37d23'30"E

then I get different result:

cct -z1000 -t0 +proj=molodensky +ellps=clrk80 +da=-112.145 +df=-0.54750714e-4 +dx=-175 +dy=-23 +dz=-303 +to +proj=latlon +datum=WGS84 file1
  5.9999485375   36.9987756806      676.6816        0.0000

what am I doing wrong there?

@kbevers
Copy link
Member

kbevers commented Sep 7, 2018

There are several things you are doing wrong.

  1. cct doesnt' understand DMS style coordinates
  2. cct expect coordinate input in the order longitude, latitude, height, time (change this with -c
  3. The +to notation is only used with cs2cs - don't use it with cct, use the command I wrote in my previous post.

You can convert DMS style coordinates to decimal degrees with cs2cs:

$ echo 6d55\'57\"S 37d23\'30\"E | cs2cs -f %.8f +proj=latlong
-6.93250000     37.39166667 0.0000000

Until cct can parse DMS coordinates you can use cs2cs as a pre-filter to cct:

echo 6d55\'57\"S 37d23\'30\"E | cs2cs -f %.8f +proj=latlong | cct -c 2,1 -t0 -z1000 +proj=molodensky +ellps=clrk80 +da=-112.145 +df=-0.54750714e-4 +dx=-175 +dy=-23 +dz=-303
 37.3924628036   -6.9346316313      991.7547        0.000

I believe we have exhausted this topic now. Please consult the documentation if you have any further problems.

@gnusupport
Copy link
Author

Thank you

@gnusupport
Copy link
Author

I am trying to tune it, as the conversion from ARC1960 is not accurate to WGS84. Please help me how to tune it so that it becomes accurate?

gie:

<gie>

operation       +proj=molodensky +ellps=clrk80 +da=-112.145 +df=-0.54750714e-4 +dx=-175 +dy=-23 +dz=-303
tolerance       2 m

accept          -1.47666 34.56861
expect          -1.47927 34.56933
# The point was accurately verified to be ARC1960 by GEOTRANS and by Tanzanian mining cadastre
</gie>

Error:

     FAILURE in gie(14):
     expected: -1.47927 34.56933
     got:      -1.476959691532   34.564311215701
     deviation:  595729.781856 mm,  expected:  2000.000000 mm
-------------------------------------------------------------------------------
total:  0 tests succeeded,  0 tests skipped,  1 tests FAILED!

how can I tune it?

@kbevers
Copy link
Member

kbevers commented Jun 1, 2019

The transformation parameters you are using are probably not the same as what is used by GEOTRANS. If you just let PROJ do it's thing you get more or less the same as you expect in your gie file:

$ echo -1.47666 34.56861 | cs2cs -f %.5f "Arc 1960" WGS84
-1.47926	34.56938 0.00000

Run the command below to find the transformation that is used. Note that several possibilites show up - choose the one that applies to your situation (likely one of the two first). You'll notice that the transformations are based on the Helmert transform and not the Molodensky as used above.

projinfo -s "Arc 1960" -t WGS84 --spatial-test intersects

@gnusupport
Copy link
Author

Thank you! I am getting almost same values as expected.

I would like to test those few options. If I just put PROJ string in gie file:
operation +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=longlat +a=6378249.145 +rf=293.465 +step +proj=push +v_3 +step +proj=cart +a=6378249.145 +rf=293.465 +step +proj=helmert +x=-160 +y=-6 +z=-302 +step +inv +pr oj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

then such is giving me error proj_create: Error -4: projection not named

How to test the PROJ strings from projinfo -s ?

@kbevers
Copy link
Member

kbevers commented Jun 2, 2019

You have a problem here: +pr oj=cart - you need to pay attention when copying and pasting.

Here's a clean version:

<gie>

operation   +proj=pipeline
            +step +proj=axisswap +order=2,1
            +step +proj=unitconvert +xy_in=deg +xy_out=rad
            +step +inv +proj=longlat +a=6378249.145 +rf=293.465
            +step +proj=push +v_3 +step +proj=cart +a=6378249.145 +rf=293.465
            +step +proj=helmert +x=-160 +y=-6 +z=-302 +step +inv
            +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step
            +proj=unitconvert +xy_in=rad +xy_out=deg +step
            +proj=axisswap +order=2,1

accept          -1.47666 34.56861
expect          -1.47927 34.56933

</gie>

@gnusupport
Copy link
Author

Oh, thank you. It looks like the last option operation +proj=noop is working closest.
Thanks much.

@kbevers
Copy link
Member

kbevers commented Jun 2, 2019

It looks like the last option operation +proj=noop is working closest.

Maybe you should take a look at the documentation for that particular operation before you do anything else ;)

https://proj.org/operations/conversions/noop.html

@gnusupport
Copy link
Author

haaaa. I get it.

I was thinking this one here is noop:
cs2cs -f %.5f "Arc 1960" WGS84

As those PROJ lines I don't know how to use. echo coord | proj projline is not how it works. Error is: can't initialize operations that take non-angular input coordinates

on

echo -1.47666 34.56861 | proj +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=longlat +a=6378249.145 +rf=293.465 +step +proj=push +v_3 +step +proj=cart +a=6378249.145 +rf=293.465 +step +proj=helmert +x=-153 +y=-5 +z=-292 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

@gnusupport
Copy link
Author

echo -1.47666 34.56861 0.0 | cct +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=longlat +a=6378249.145 +rf=293.465 +step +proj=push +v_3 +step +proj=cart +a=6378249.145 +rf=293.465 +step +proj=helmert +x=-153 +y=-5 +z=-292 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

This works, so this way I can test lines.

@kbevers
Copy link
Member

kbevers commented Jun 2, 2019

Generally cs2cs -f %.5f "Arc 1960" WGS84 should give you the correct transformation. It takes your input data into consideration and chooses the best transformation based on the area of use for the different transformations between the two reference systems. I would use that if I were you.

Note that your input has to be served as latitude, longitude and not the opposite way round. This is because cs2cs respects the axis order that the local geodetic authority has defined. When in doubt you can always check the output of projinfo "Arc 1960" to see how the axes are defined. If you inspect the pipelines above you will notice that the first and last step changes the order of the axis.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants