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.cos;
|
---|
9 | import static java.lang.Math.exp;
|
---|
10 | import static java.lang.Math.log;
|
---|
11 | import static java.lang.Math.pow;
|
---|
12 | import static java.lang.Math.sin;
|
---|
13 | import static java.lang.Math.sqrt;
|
---|
14 | import static java.lang.Math.tan;
|
---|
15 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
16 |
|
---|
17 | import org.openstreetmap.josm.data.Bounds;
|
---|
18 | import org.openstreetmap.josm.data.projection.Ellipsoid;
|
---|
19 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
20 | import org.openstreetmap.josm.tools.Utils;
|
---|
21 |
|
---|
22 | /**
|
---|
23 | * Implementation of the stereographic double projection,
|
---|
24 | * also known as Oblique Stereographic and the Schreiber double stereographic projection.
|
---|
25 | *
|
---|
26 | * @author vholten
|
---|
27 | *
|
---|
28 | * Source: IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2,
|
---|
29 | * Sec. 1.3.7.1 Oblique and Equatorial Stereographic, http://www.epsg.org/GuidanceNotes
|
---|
30 | */
|
---|
31 | public class DoubleStereographic extends AbstractProj {
|
---|
32 |
|
---|
33 | private Ellipsoid ellps;
|
---|
34 | private double n;
|
---|
35 | private double c;
|
---|
36 | private double chi0;
|
---|
37 | private double r;
|
---|
38 |
|
---|
39 | private static final double EPSILON = 1e-12;
|
---|
40 |
|
---|
41 | @Override
|
---|
42 | public String getName() {
|
---|
43 | return tr("Double Stereographic");
|
---|
44 | }
|
---|
45 |
|
---|
46 | @Override
|
---|
47 | public String getProj4Id() {
|
---|
48 | return "sterea";
|
---|
49 | }
|
---|
50 |
|
---|
51 | @Override
|
---|
52 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
53 | super.initialize(params);
|
---|
54 | if (params.lat0 == null)
|
---|
55 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
|
---|
56 | ellps = params.ellps;
|
---|
57 | initialize(params.lat0);
|
---|
58 | }
|
---|
59 |
|
---|
60 | private void initialize(double lat0) {
|
---|
61 | double phi0 = Utils.toRadians(lat0);
|
---|
62 | double e2 = ellps.e2;
|
---|
63 | r = sqrt(1-e2) / (1 - e2*pow(sin(phi0), 2));
|
---|
64 | n = sqrt(1 + ellps.eb2 * pow(cos(phi0), 4));
|
---|
65 | double s1 = (1 + sin(phi0)) / (1 - sin(phi0));
|
---|
66 | double s2 = (1 - e * sin(phi0)) / (1 + e * sin(phi0));
|
---|
67 | double w1 = pow(s1 * pow(s2, e), n);
|
---|
68 | double sinchi00 = (w1 - 1) / (w1 + 1);
|
---|
69 | c = (n + sin(phi0)) * (1 - sinchi00) / ((n - sin(phi0)) * (1 + sinchi00));
|
---|
70 | double w2 = c * w1;
|
---|
71 | chi0 = asin((w2 - 1) / (w2 + 1));
|
---|
72 | }
|
---|
73 |
|
---|
74 | @Override
|
---|
75 | public double[] project(double phi, double lambda) {
|
---|
76 | double nLambda = n * lambda;
|
---|
77 | double sa = (1 + sin(phi)) / (1 - sin(phi));
|
---|
78 | double sb = (1 - e * sin(phi)) / (1 + e * sin(phi));
|
---|
79 | double w = c * pow(sa * pow(sb, e), n);
|
---|
80 | double chi = asin((w - 1) / (w + 1));
|
---|
81 | double b = 1 + sin(chi) * sin(chi0) + cos(chi) * cos(chi0) * cos(nLambda);
|
---|
82 | double x = 2 * r * cos(chi) * sin(nLambda) / b;
|
---|
83 | double y = 2 * r * (sin(chi) * cos(chi0) - cos(chi) * sin(chi0) * cos(nLambda)) / b;
|
---|
84 | return new double[] {x, y};
|
---|
85 | }
|
---|
86 |
|
---|
87 | @Override
|
---|
88 | public double[] invproject(double x, double y) {
|
---|
89 | double e2 = ellps.e2;
|
---|
90 | double g = 2 * r * tan(PI/4 - chi0/2);
|
---|
91 | double h = 4 * r * tan(chi0) + g;
|
---|
92 | double i = atan(x/(h + y));
|
---|
93 | double j = atan(x/(g - y)) - i;
|
---|
94 | double chi = chi0 + 2 * atan((y - x * tan(j/2)) / (2 * r));
|
---|
95 | double lambda = (j + 2*i) / n;
|
---|
96 | double psi = 0.5 * log((1 + sin(chi)) / (c*(1 - sin(chi)))) / n;
|
---|
97 | double phiprev = -1000;
|
---|
98 | int iteration = 0;
|
---|
99 | double phi = 2 * atan(exp(psi)) - PI/2;
|
---|
100 | while (abs(phi - phiprev) > EPSILON) {
|
---|
101 | if (++iteration > 10)
|
---|
102 | throw new IllegalStateException("Too many iterations");
|
---|
103 | phiprev = phi;
|
---|
104 | double psii = log(tan(phi/2 + PI/4) * pow((1 - e * sin(phi)) / (1 + e * sin(phi)), e/2));
|
---|
105 | phi = phi - (psii - psi) * cos(phi) * (1 - e2 * pow(sin(phi), 2)) / (1 - e2);
|
---|
106 | }
|
---|
107 | return new double[] {phi, lambda};
|
---|
108 | }
|
---|
109 |
|
---|
110 | @Override
|
---|
111 | public Bounds getAlgorithmBounds() {
|
---|
112 | return new Bounds(-89, -87, 89, 87, false);
|
---|
113 | }
|
---|
114 | }
|
---|