1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import static java.lang.Math.PI;
|
---|
5 | import static java.lang.Math.abs;
|
---|
6 | import static java.lang.Math.asin;
|
---|
7 | import static java.lang.Math.atan;
|
---|
8 | import static java.lang.Math.atan2;
|
---|
9 | import static java.lang.Math.cos;
|
---|
10 | import static java.lang.Math.exp;
|
---|
11 | import static java.lang.Math.log;
|
---|
12 | import static java.lang.Math.pow;
|
---|
13 | import static java.lang.Math.sin;
|
---|
14 | import static java.lang.Math.sqrt;
|
---|
15 | import static java.lang.Math.tan;
|
---|
16 | import static java.lang.Math.toRadians;
|
---|
17 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
18 |
|
---|
19 | import org.openstreetmap.josm.data.projection.Ellipsoid;
|
---|
20 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
21 |
|
---|
22 | // CHECKSTYLE.OFF: LineLength
|
---|
23 |
|
---|
24 | /**
|
---|
25 | * Projection for the SwissGrid CH1903 / L03, see <a href="https://en.wikipedia.org/wiki/Swiss_coordinate_system">Wikipedia article</a>.<br>
|
---|
26 | *
|
---|
27 | * Calculations were originally based on
|
---|
28 | * <a href="http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.12749.DownloadFile.tmp/ch1903wgs84en.pdf">
|
---|
29 | * simple formula</a>.<br>
|
---|
30 | *
|
---|
31 | * August 2010 update to
|
---|
32 | * <a href="http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf">
|
---|
33 | * this formula (rigorous formulas)</a>.
|
---|
34 | */
|
---|
35 | public class SwissObliqueMercator implements Proj {
|
---|
36 |
|
---|
37 | // CHECKSTYLE.ON: LineLength
|
---|
38 |
|
---|
39 | private Ellipsoid ellps;
|
---|
40 | private double kR;
|
---|
41 | private double alpha;
|
---|
42 | private double b0;
|
---|
43 | private double k;
|
---|
44 |
|
---|
45 | private static final double EPSILON = 1e-11;
|
---|
46 |
|
---|
47 | @Override
|
---|
48 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
49 | if (params.lat0 == null)
|
---|
50 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
|
---|
51 | ellps = params.ellps;
|
---|
52 | initialize(params.lat0);
|
---|
53 | }
|
---|
54 |
|
---|
55 | private void initialize(double lat_0) {
|
---|
56 | double phi0 = toRadians(lat_0);
|
---|
57 | kR = sqrt(1 - ellps.e2) / (1 - (ellps.e2 * pow(sin(phi0), 2)));
|
---|
58 | alpha = sqrt(1 + (ellps.eb2 * pow(cos(phi0), 4)));
|
---|
59 | b0 = asin(sin(phi0) / alpha);
|
---|
60 | k = log(tan(PI / 4 + b0 / 2)) - alpha
|
---|
61 | * log(tan(PI / 4 + phi0 / 2)) + alpha * ellps.e / 2
|
---|
62 | * log((1 + ellps.e * sin(phi0)) / (1 - ellps.e * sin(phi0)));
|
---|
63 | }
|
---|
64 |
|
---|
65 | @Override
|
---|
66 | public String getName() {
|
---|
67 | return tr("Swiss Oblique Mercator");
|
---|
68 | }
|
---|
69 |
|
---|
70 | @Override
|
---|
71 | public String getProj4Id() {
|
---|
72 | return "somerc";
|
---|
73 | }
|
---|
74 |
|
---|
75 | @Override
|
---|
76 | public double[] project(double phi, double lambda) {
|
---|
77 |
|
---|
78 | double S = alpha * log(tan(PI / 4 + phi / 2)) - alpha * ellps.e / 2
|
---|
79 | * log((1 + ellps.e * sin(phi)) / (1 - ellps.e * sin(phi))) + k;
|
---|
80 | double b = 2 * (atan(exp(S)) - PI / 4);
|
---|
81 | double l = alpha * lambda;
|
---|
82 |
|
---|
83 | double lb = atan2(sin(l), sin(b0) * tan(b) + cos(b0) * cos(l));
|
---|
84 | double bb = asin(cos(b0) * sin(b) - sin(b0) * cos(b) * cos(l));
|
---|
85 |
|
---|
86 | double y = kR * lb;
|
---|
87 | double x = kR / 2 * log((1 + sin(bb)) / (1 - sin(bb)));
|
---|
88 |
|
---|
89 | return new double[] {y, x};
|
---|
90 | }
|
---|
91 |
|
---|
92 | @Override
|
---|
93 | public double[] invproject(double y, double x) {
|
---|
94 | double lb = y / kR;
|
---|
95 | double bb = 2 * (atan(exp(x / kR)) - PI / 4);
|
---|
96 |
|
---|
97 | double b = asin(cos(b0) * sin(bb) + sin(b0) * cos(bb) * cos(lb));
|
---|
98 | double l = atan2(sin(lb), cos(b0) * cos(lb) - sin(b0) * tan(bb));
|
---|
99 |
|
---|
100 | double lambda = l / alpha;
|
---|
101 | double phi = b;
|
---|
102 | double s = 0;
|
---|
103 |
|
---|
104 | double prevPhi = -1000;
|
---|
105 | int iteration = 0;
|
---|
106 | // iteration to finds S and phi
|
---|
107 | while (abs(phi - prevPhi) > EPSILON) {
|
---|
108 | if (++iteration > 30)
|
---|
109 | throw new RuntimeException("Two many iterations");
|
---|
110 | prevPhi = phi;
|
---|
111 | s = 1 / alpha * (log(tan(PI / 4 + b / 2)) - k) + ellps.e
|
---|
112 | * log(tan(PI / 4 + asin(ellps.e * sin(phi)) / 2));
|
---|
113 | phi = 2 * atan(exp(s)) - PI / 2;
|
---|
114 | }
|
---|
115 | return new double[] {phi, lambda};
|
---|
116 | }
|
---|
117 | }
|
---|