Opened 3 years ago

Closed 3 years ago

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

Reported by: Owned by: sbrunner team normal Core patch

### Description

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

### Changed 3 years ago by sbrunner

Patch to use the rigorous formulas

### comment:1 follow-up: ↓ 2 Changed 3 years ago by bastiK

Have you tested against some coordinates given in that document?

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

### comment:2 in reply to: ↑ 1 Changed 3 years ago by anonymous

Hi,

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

### comment:3 Changed 3 years ago by 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

### comment:4 Changed 3 years ago by anonymous

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

CU
Stéphane

Patch with atan2

### comment:5 Changed 3 years ago by 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

### comment:6 Changed 3 years ago by bastiK

• Summary changed from [patch] Swissgrild is unusable because he use aprximate formulas => rigorous formulas to [patch] Swissgrild uses approximate formulas => better use exact formulas

### comment:7 Changed 3 years ago by 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.

### comment:8 Changed 3 years ago by bastiK

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

### comment:9 Changed 3 years ago by 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

### comment:10 Changed 3 years ago by 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.

### comment:11 Changed 3 years ago by 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.

### comment:12 Changed 3 years ago by 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. :)

### comment:13 Changed 3 years ago by bastiK

• 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?

### comment:14 follow-up: ↓ 15 Changed 3 years ago by 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

### Changed 3 years ago by sbrunner

Manly use Ellipsoide.

### comment:15 in reply to: ↑ 14 ; follow-up: ↓ 16 Changed 3 years ago by bastiK

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 !

• 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.

### comment:16 in reply to: ↑ 15 ; follow-up: ↓ 17 Changed 3 years ago by bastiK

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 !

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.

### comment:17 in reply to: ↑ 16 ; follow-up: ↓ 18 Changed 3 years ago by 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 :

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

### Changed 3 years ago by sbrunner

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

### comment:18 in reply to: ↑ 17 Changed 3 years ago by 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 :

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

### comment:19 Changed 3 years ago by bastiK

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

### comment:20 follow-up: ↓ 21 Changed 3 years ago by 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+.

### comment:21 in reply to: ↑ 20 ; follow-up: ↓ 22 Changed 3 years ago by sbrunner

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 :

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

### comment:22 in reply to: ↑ 21 ; follow-up: ↓ 23 Changed 3 years ago by 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 :

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...

### comment:23 in reply to: ↑ 22 ; follow-up: ↓ 24 Changed 3 years ago by 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

### comment:24 in reply to: ↑ 23 ; follow-up: ↓ 26 Changed 3 years ago by anonymous

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.)

anonymous <- me

### comment:26 in reply to: ↑ 24 Changed 3 years ago by 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.

### comment:27 Changed 3 years ago by 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.

### Changed 3 years ago by sbrunner

decrase pressision requirement, a fix

### comment:28 Changed 3 years ago by bastiK

• Type changed from defect to enhancement

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.

### comment:30 Changed 3 years ago by 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.

### comment:31 Changed 3 years ago by bastiK

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

### comment:32 Changed 3 years ago by 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?

### comment:33 Changed 3 years ago by jttt

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

### comment:34 in reply to: ↑ 29 Changed 3 years ago by anonymous

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 !

### comment:35 Changed 3 years ago by sbrunner

• Resolution set to fixed
• Status changed from new to closed

anonymous in me !

### comment:36 Changed 3 years ago by bastiK

Does that help?

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

### comment:37 Changed 3 years ago by 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.

### comment:38 Changed 3 years ago by bastiK

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

### Modify Ticket

Change Properties