[patch] Swissgrild uses approximate formulas => better use exact formulas

[sbrunner]

In the swisstopo there is mentioned two methods to calculate the swissgrid, actually it's the approximate method who is implemented.

​http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf

[sbrunner]

Patch to use the rigorous formulas

[bastiK]

Have you tested against some coordinates given in that document?

(Maybe it is saver to use atan2 instead of atan?)

[anonymous]

Hi,

Replying to bastiK:

Have you tested against some coordinates given in that document?

Yes I wave testes with some coordinates in swiss extremities, and compare the result with the new method and the old one (to escape to be worth !)

And finally I test it with the epfl wms :
​http://wiki.openstreetmap.org/wiki/EPFL_WMS.

(Maybe it is saver to use atan2 instead of atan?)

I don't know but I don't understand what the atan2 do ... ?

CU
Stéphane

[stoecker]

Basically atan2(x,y) is atan(x/y), but cares for the correct sector when x and/or y are negative: ​http://en.wikipedia.org/wiki/Atan2

[anonymous]

Ok, thanks, Than I will test it this evening ;-)

CU
Stéphane

[sbrunner]

Patch with atan2

[sbrunner]

Hello,

I just use the atan2 method, and effectively when I see the formulas who use the atan it seem relay useful ;-)

CU
Stéphane

[bastiK]

Not so good to do a fixed of number of iterations. Better add a proper terminating condition and limit the maximum number of iterations. Throw exception if it hasn't terminated by then.

Also, please clean up the patch.

[bastiK]

If you like, you can add test cases. (See folder test/ for examples.)

[sbrunner]

cleaning, better iterantion, add test

[sbrunner]

Hello,

I just change the iteration as you say, move the test as unit test and put only the needed file in the patch (my previous patch is really poor :( )

CU
Stéphane

[bastiK]

Ok, I don't want to be too nitpicking, because I have been already :), but 2 things:

If I remember correctly, there are 2 variables in the iteration, so it should be

while (Math.abs(phi - prevPhi) > 0.0000001 || Math.abs(Rn(?) - prevRn) > EPSILON) {

I would insert it myself, but I'm not sure, what is a reasonable value for EPSILON?

If the result covers a larger range, I would write abs(phi - prevPhi) > EPS * max(abs(phi),abs(prevPhi)) or something like that.

E.g. the current condition is fulfilled for phi=0.0 & prevPhi=0.00000009 but not for
phi=1000000.0000001 & prevPhi=phi=1000000.0000003.

[bastiK]

Last part should be:

If the result covers a larger range, I would write abs(phi - prevPhi) > EPS * min(abs(phi),abs(prevPhi)) or something like that.

E.g. the current condition is fulfilled for phi=0.0 & prevPhi=0.00000009 but not for phi=1000000.0000001 & prevPhi=1000000.0000003.

[bastiK]

Damn, what I really mean is this part:

  // 5 iteration to fins S and phi 
 	        for (int i = 0; i < 6; i++) { 
 	            S = 1 / alpha * (Math.log(Math.tan(Math.PI / 4 + b / 2)) - K) + E 
 	                    * Math.log(Math.tan(Math.PI / 4 + Math.asin(E * Math.sin(phi)) / 2)); 
 	            phi = 2 * Math.atan(Math.exp(S)) - Math.PI / 2; 
 	        } 

Both S and phi must change only little at the same time. Better not guess that 5 iterations are enough...

I would like to look at it in more detail, but not today. :)

[bastiK]

More comments:

• There is a method cart2LatLon in UTM_France_DOM.java; It should be used if possible (and maybe moved to Ellipoide.java).
• I know WSG 84 and GRS 80 are almost identical, but why mix them up unnecessarily? (It is called correctETRS89toCH1903(), then you use WSG constants.
• Why do you discard the height information when translating from WSG to Bessel?

[anonymous]

Hello,

Many things ;-)

• For the phi and Rn I don't thing tat it useful to to test a temporary variable if th value of phi don't change it must mean that the value of phi don't change any more.

• For the test of the delta phi I don't agree with you because phi is an angle than we shouldn't be more precise near the 0 value than near th pi value !

• For the 5 iteration to fins S and phi => oups I miss it I probably to tied last evening ! => corrected

• I move the cart2LatLon method to Ellipoide an add the revere one. I use the ant its relay clearer !

• The Ellipoide use by the swiss grid is real different If I miss it I have some error of just less than 100m !!! I thing that with the code that use the Ellipoide class it's clearer.

• For the height it's just because OSM and JOSM don't use it ! (is it a mistake from me ?)

CU
Stéphane

[sbrunner]

Manly use Ellipsoide.

[bastiK]

Replying to anonymous:

Hello,

Many things ;-)

• For the phi and Rn I don't thing tat it useful to to test a temporary variable if th value of phi don't change it must mean that the value of phi don't change any more.

Yes, I mixed it up with phi & S

• For the test of the delta phi I don't agree with you because phi is an angle than we shouldn't be more precise near the 0 value than near th pi value !

True, what about S?

• For the 5 iteration to fins S and phi => oups I miss it I probably to tied last evening ! => corrected

• I move the cart2LatLon method to Ellipoide an add the revere one. I use the ant its relay clearer !

For consistency, LatLon in JOSM should be in deg and not rad. (I know, it is wrong in UTM_France_DOM.java, as well.)

• The Ellipoide use by the swiss grid is real different If I miss it I have some error of just less than 100m !!! I thing that with the code that use the Ellipoide class it's clearer.

What you call SWISS_ELLIPSOID is actually Bessel ellipsoid. Also, why do you write ETRS89? This system has nothing to do with the conversion CH1903+ <-> WSG 84, right?

There seems to be a difference between CH1903 and CH1903+. Can you clarify in the comments?

• For the height it's just because OSM and JOSM don't use it ! (is it a mistake from me ?)

Right, I thougt h would enter in the formula for swiss coordinates.

[bastiK]

Replying to bastiK:

Replying to anonymous:

Hello,

Many things ;-)

• For the phi and Rn I don't thing tat it useful to to test a temporary variable if th value of phi don't change it must mean that the value of phi don't change any more.

Yes, I mixed it up with phi & S

• For the test of the delta phi I don't agree with you because phi is an angle than we shouldn't be more precise near the 0 value than near th pi value !

True, what about S?

Forget what I said, S is also just a auxiliary parameter, so there is nothing wrong here, sorry.

There seems to be a difference between CH1903 and CH1903+. Can you clarify in the comments?

I looked it up: It seems that we still not have sub-meter accuracy:

• What you are calculating is CH1903+ but without the offset of 2000000 & 1000000
• To get CH1903 (which seems to be still in use), you need triangulation correction (so called FINELTRA method). The corrections are up to 3 meters.
• I cannot download the CHENyx06 file from their site, but I'm quite sure it is too large to be distributed along with JOSM core. The simplified NTv2 grid is ~1MB (compressed). However, a plugin would be fine.

[sbrunner]

What you call SWISS_ELLIPSOID is actually Bessel ellipsoid. Also, why do you write ETRS89? This system has nothing to do with the conversion CH1903+ <-> WSG 84, right?

Effectively I rename them.

There seems to be a difference between CH1903 and CH1903+. Can you clarify in the comments?

As I understand the CH1903+ it just small improvement but I don't know exactly what !

For the height it's just because OSM and JOSM don't use it ! (is it a mistake from me ?)

Right, I thougt h would enter in the formula for swiss coordinates.

Effectively it's because the height is not the same in CH1903+ than in WSG84 in the example there is 5 cm of difference.

I looked it up: It seems that we still not have sub-meter accuracy:

In my test I have a precision of about 5cm than it's more precise than OSM !

And the big deference than before is that the booth transformation send and return has an error that's is less than 2mm !

• What you are calculating is CH1903+ but without the offset of 2000000 & 1000000

Effectively I use the LV03 coordinates, not the LV95 as all service I use, for examples :

• ​http://map.geo.admin.ch
• ​http://plan.epfl.ch
• ​http://wiki.openstreetmap.org/wiki/EPFL_WMS

And the previous version will do the same.

• To get CH1903 (which seems to be still in use), you need triangulation correction (so called FINELTRA method). The corrections are up to 3 meters.
• I cannot download the CHENyx06 file from their site, but I'm quite sure it is too large to be distributed along with JOSM core. The simplified NTv2 grid is ~1MB (compressed). However, a plugin would be fine.

It's possible but I don't think that we need it.

CU
Stéphane

[sbrunner]

Rename come things, use degree in methode of Ellipsoid object.

[sbrunner]

• What you are calculating is CH1903+ but without the offset of 2000000 & 1000000

Effectively I use the LV03 coordinates, not the LV95 as all service I use, for examples :

• ​http://map.geo.admin.ch
• ​http://plan.epfl.ch
• ​http://wiki.openstreetmap.org/wiki/EPFL_WMS

And the previous version will do the same.

The LV95 is the projection has EPSG:2056
​http://spatialreference.org/ref/epsg/2056/

CU
Stéphane

[bastiK]

(In [3473]) see #5327 (patch by sbrunner) - Swissgrild uses approximate formulas => better use exact formulas

[bastiK]

Where is your test data from? I added the examples from the swisstopo pdf and it demonstrates the local distortions of CH1903 compared to CH1903+.

[sbrunner]

Replying to bastiK:

Where is your test data from? I added the examples from the swisstopo pdf and it demonstrates the local distortions of CH1903 compared to CH1903+.

I base on the exact result specified on page 11 and 12 of the reference document :
​http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf.

Thanks for the checkin, I see that you increase the precision of phi (DELTA_PHI) from 1e-8 to 1e-11, I don't know that its useful because when I don some test on send and return calculation the prescription don't increase with value lower than 1e-8. I say that but finally it's not important for me ;-)

CU
Stéphane

[bastiK]

Replying to sbrunner:

Replying to bastiK:

Where is your test data from? I added the examples from the swisstopo pdf and it demonstrates the local distortions of CH1903 compared to CH1903+.

I base on the exact result specified on page 11 and 12 of the reference document :
​http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf.

When I open this document, on page 11 and 12 there is no such information. I'm looking for the source of the coordinates for Lausanne, Schaffouse,...

However on p. 17 there is data about Zimmerwald, Chrischona, ... which I included in the test file. These tests fail at the moment, but it seems, that is expected and nothing to worry about.

Thanks for the checkin, I see that you increase the precision of phi (DELTA_PHI) from 1e-8 to 1e-11, I don't know that its useful because when I don some test on send and return calculation the prescription don't increase with value lower than 1e-8. I say that but finally it's not important for me ;-)

I did it, because 1e-11 was already used in other places, but it probably does not matter...

[sbrunner]

When I open this document, on page 11 and 12 there is no such information. I'm looking for the source of the coordinates for Lausanne, Schaffouse,...

Effectively only the Ref will be taken on this document I put the other to prevent side effect when I change something in the patch and now when somebody change something in the ellipsoid.

However on p. 17 there is data about Zimmerwald, Chrischona, ... which I included in the test file. These tests fail at the moment, but it seems, that is expected and nothing to worry about.

On page 17 they don't specify the WSG84 coordinate than it can be directly used do do a test..

I did it, because 1e-11 was already used in other places, but it probably does not matter...

OK, no problem for me ;-)

CU
Stéphane

[anonymous]

Replying to sbrunner:

When I open this document, on page 11 and 12 there is no such information. I'm looking for the source of the coordinates for Lausanne, Schaffouse,...

Effectively only the Ref will be taken on this document I put the other to prevent side effect when I change something in the patch and now when somebody change something in the ellipsoid.

Well, it is not a really good test for a calculation to compare with the results of that very same calculation. ;)

At least it should be noted, that it is not independent data.

However on p. 17 there is data about Zimmerwald, Chrischona, ... which I included in the test file. These tests fail at the moment, but it seems, that is expected and nothing to worry about.

On page 17 they don't specify the WSG84 coordinate than it can be directly used do do a test..

GRS 80 it is practically the same as WSG 84. The differences in the projected coordinates are in the nanometer scale. (The distinction is only relevant for astronomers.)

[bastiK]

anonymous <- me

[sbrunner]

However on p. 17 there is data about Zimmerwald, Chrischona, ... which I included in the test file. These tests fail at the moment, but it seems, that is expected and nothing to worry about.

On page 17 they don't specify the WSG84 coordinate than it can be directly used do do a test..

GRS 80 it is practically the same as WSG 84. The differences in the projected coordinates are in the nanometer scale. (The distinction is only relevant for astronomers.)

It seem that your true, I will have a look on this this evening.

[sbrunner]

Hello,

Effectively it seem that on the reference purpose on the official document we need the chenyx06 data set to have sub-metric precision.

Finally the real advantage than before is when we do a conversion in booth side (send and return) we have nearly the same result.

CU
Stéphane

PS: it you want to make the test passing I purpose a patch who decrease the precession requirement, and do a fix on missing parenthesis.

[sbrunner]

decrase pressision requirement, a fix

[bastiK]

Hi, 1-2 meters is more or less what I would consider the threshold for being accurate enough for osm use. If the tests fail, that is fine, I can live with that. Maybe it motivates someone to get the chenyx06.bin file and fix it. :)

Also would be nice to add CH1903+ support as it seems to be the upcoming standard.

Anyway, feel free to close the ticket, if you think it is Ok like this.

[bastiK]

You can download the CHENYX06.gsb file here:

​http://trac.osgeo.org/csmap/browser/trunk/CsMapDev/Dictionaries/Switzerland

[jttt]

But I can't live with failing tests :-). I'm running hudson on my laptop and I get e-mail everytime some test fails. I might not notice some real problems when projection tests will always fail.

[bastiK]

(In [3483]) fixed tests for projection (see #5327)

[bastiK]

I added another unrelated test that happens to fail at the moment.

There is a System property to suppress these errors, now. Is that Ok for you?

[jttt]

Well, I usually simply use @Ignore annotation on tests that are known to be broken, but system property is also ok.

[anonymous]

Replying to bastiK:

You can download the CHENYX06.gsb file here:

​http://trac.osgeo.org/csmap/browser/trunk/CsMapDev/Dictionaries/Switzerland

Effectively, we can also download it from here - ​http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/lv03-lv95/chenyx06/distortion_grids.html but I don't know how to use it, the notice in only in German and I speak only French and English !

[sbrunner]

anonymous in me !

[bastiK]

Does that help?

​http://www.swisstopo.admin.ch/internet/swisstopo/en/home/products/software/products/chenyx06.parsysrelated1.72268.downloadList.91729.DownloadFile.tmp/fineltrafr.pdf

If you need a German document translated and google does not help, I'll see what I can do.

[bastiK]

All files I have downloaded so far (CHENyx06.gsb, CHENyx06a.gsb) are grid interpolations in standard NTv2 format. Note, that JOSM already supports this format, see NTV2GridShift.java.

But it would be nice to have the real triangulation data, which is described in the document I linked above.

[bastiK]

Ok, let's finish here, I opened a new ticket for the remaining things: #5399