source: josm/trunk/src/org/openstreetmap/josm/data/projection/proj/SwissObliqueMercator.java@ 4382

Last change on this file since 4382 was 4285, checked in by bastiK, 13 years ago

major projection rework

More modular structure, inspired by Proj.4.

There are almost no semantic changes to the projection algorithms. Mostly factors of 'a' and 180/PI have been moved from one place to the other. In UTM_France_DOM, the ellipsoid conversion for the first 3 projections has been changed from hayford <-> GRS80 to hayford <-> WGS84.

Some redundant algorithms have been removed. In particular:

  • UTM_France_DOM used to have its own Transverse Mercator implementation. It is different from the implementation in TransverseMercator.java as it has another series expansion. For EPSG::2975, there are numeric differences on centimeter scale. However, the new data fits better with Proj.4 output.
  • Also removed are alternate implementations of LambertConformalConic. (They are all quite similar, though.)
  • Property svn:eol-style set to native
File size: 3.2 KB
Line 
1//License: GPL. For details, see LICENSE file.
2package org.openstreetmap.josm.data.projection.proj;
3
4import static java.lang.Math.*;
5
6import static org.openstreetmap.josm.tools.I18n.tr;
7
8import org.openstreetmap.josm.data.projection.Ellipsoid;
9
10/**
11 * Projection for the SwissGrid CH1903 / L03, see http://de.wikipedia.org/wiki/Swiss_Grid.
12 *
13 * Calculations are based on formula from
14 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf
15 *
16 * August 2010 update to this formula (rigorous formulas)
17 * http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf
18 */
19public class SwissObliqueMercator implements Proj {
20
21 private final Ellipsoid ellps;
22 private double kR;
23 private double alpha;
24 private double b0;
25 private double K;
26
27 private static final double EPSILON = 1e-11;
28
29 public SwissObliqueMercator(Ellipsoid ellps, double lat_0) {
30 this.ellps = ellps;
31 updateParameters(lat_0);
32 }
33
34 public void updateParameters(double lat_0) {
35 double phi0 = toRadians(lat_0);
36 kR = sqrt(1 - ellps.e2) / (1 - (ellps.e2 * pow(sin(phi0), 2)));
37 alpha = sqrt(1 + (ellps.eb2 * pow(cos(phi0), 4)));
38 b0 = asin(sin(phi0) / alpha);
39 K = log(tan(PI / 4 + b0 / 2)) - alpha
40 * log(tan(PI / 4 + phi0 / 2)) + alpha * ellps.e / 2
41 * log((1 + ellps.e * sin(phi0)) / (1 - ellps.e * sin(phi0)));
42 }
43
44 @Override
45 public String getName() {
46 return tr("Swiss Oblique Mercator");
47 }
48
49 @Override
50 public String getProj4Id() {
51 return "somerc";
52 }
53
54 @Override
55 public double[] project(double phi, double lambda) {
56
57 double S = alpha * log(tan(PI / 4 + phi / 2)) - alpha * ellps.e / 2
58 * log((1 + ellps.e * sin(phi)) / (1 - ellps.e * sin(phi))) + K;
59 double b = 2 * (atan(exp(S)) - PI / 4);
60 double l = alpha * lambda;
61
62 double lb = atan2(sin(l), sin(b0) * tan(b) + cos(b0) * cos(l));
63 double bb = asin(cos(b0) * sin(b) - sin(b0) * cos(b) * cos(l));
64
65 double y = kR * lb;
66 double x = kR / 2 * log((1 + sin(bb)) / (1 - sin(bb)));
67
68 return new double[] { y, x };
69 }
70
71 @Override
72 public double[] invproject(double y, double x) {
73 double lb = y / kR;
74 double bb = 2 * (atan(exp(x / kR)) - PI / 4);
75
76 double b = asin(cos(b0) * sin(bb) + sin(b0) * cos(bb) * cos(lb));
77 double l = atan2(sin(lb), cos(b0) * cos(lb) - sin(b0) * tan(bb));
78
79 double lambda = l / alpha;
80 double phi = b;
81 double S = 0;
82
83 double prevPhi = -1000;
84 int iteration = 0;
85 // iteration to finds S and phi
86 while (abs(phi - prevPhi) > EPSILON) {
87 if (++iteration > 30)
88 throw new RuntimeException("Two many iterations");
89 prevPhi = phi;
90 S = 1 / alpha * (log(tan(PI / 4 + b / 2)) - K) + ellps.e
91 * log(tan(PI / 4 + asin(ellps.e * sin(phi)) / 2));
92 phi = 2 * atan(exp(S)) - PI / 2;
93 }
94 return new double[] { phi, lambda };
95 }
96
97}
Note: See TracBrowser for help on using the repository browser.