Modify

Opened 11 years ago

Closed 11 years ago

Last modified 11 years ago

#5327 closed enhancement (fixed)

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

Reported by: sbrunner Owned by: team
Priority: normal Milestone:
Component: Core Version:
Keywords: patch Cc:

Description

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

Attachments (6)

swiss.patch (8.8 KB) - added by sbrunner 11 years ago.
Patch to use the rigorous formulas
swissgrid.patch (40.7 KB) - added by sbrunner 11 years ago.
Patch with atan2
new.patch (15.0 KB) - added by sbrunner 11 years ago.
cleaning, better iterantion, add test
swissgrid.2.patch (17.7 KB) - added by sbrunner 11 years ago.
Manly use Ellipsoide.
swissgrid.3.patch (17.9 KB) - added by sbrunner 11 years ago.
Rename come things, use degree in methode of Ellipsoid object.
swissgridtest.patch (3.6 KB) - added by sbrunner 11 years ago.
decrase pressision requirement, a fix

Download all attachments as: .zip

Change History (44)

Changed 11 years ago by sbrunner

Attachment: swiss.patch added

Patch to use the rigorous formulas

comment:1 Changed 11 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 11 years ago by 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

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

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

CU
Stéphane

Changed 11 years ago by sbrunner

Attachment: swissgrid.patch added

Patch with atan2

comment:5 Changed 11 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 11 years ago by bastiK

Summary: [patch] Swissgrild is unusable because he use aprximate formulas => rigorous formulas[patch] Swissgrild uses approximate formulas => better use exact formulas

comment:7 Changed 11 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 11 years ago by bastiK

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

Changed 11 years ago by sbrunner

Attachment: new.patch added

cleaning, better iterantion, add test

comment:9 Changed 11 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 11 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 11 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 11 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 11 years ago by 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?

comment:14 Changed 11 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 11 years ago by sbrunner

Attachment: swissgrid.2.patch added

Manly use Ellipsoide.

comment:15 in reply to:  14 ; Changed 11 years ago by 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.

comment:16 in reply to:  15 ; Changed 11 years ago by 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.

comment:17 in reply to:  16 ; Changed 11 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 11 years ago by sbrunner

Attachment: swissgrid.3.patch added

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

comment:18 in reply to:  17 Changed 11 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 11 years ago by bastiK

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

comment:20 Changed 11 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 ; Changed 11 years ago by 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

comment:22 in reply to:  21 ; Changed 11 years ago by 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...

comment:23 in reply to:  22 ; Changed 11 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 ; Changed 11 years ago by 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.)

comment:25 Changed 11 years ago by bastiK

anonymous <- me

comment:26 in reply to:  24 Changed 11 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 11 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 11 years ago by sbrunner

Attachment: swissgridtest.patch added

decrase pressision requirement, a fix

comment:28 Changed 11 years ago by bastiK

Type: defectenhancement

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:29 Changed 11 years ago by bastiK

comment:30 Changed 11 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 11 years ago by bastiK

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

comment:32 Changed 11 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 11 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 11 years ago by 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 !

comment:35 Changed 11 years ago by sbrunner

Resolution: fixed
Status: newclosed

anonymous in me !

comment:36 Changed 11 years ago by bastiK

comment:37 Changed 11 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 11 years ago by bastiK

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

Modify Ticket

Change Properties
Set your email in Preferences
Action
as closed The owner will remain team.
as The resolution will be set.
The resolution will be deleted.

Add Comment


E-mail address and name can be saved in the Preferences.

 
Note: See TracTickets for help on using tickets.