1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
5 |
|
---|
6 | import org.openstreetmap.josm.data.Bounds;
|
---|
7 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
8 |
|
---|
9 | /**
|
---|
10 | * The polar case of the stereographic projection.
|
---|
11 | * <p>
|
---|
12 | * In the proj.4 library, the code "stere" covers several variants of the
|
---|
13 | * Stereographic projection, depending on the latitude of natural origin
|
---|
14 | * (parameter lat_0).
|
---|
15 | * <p>
|
---|
16 | *
|
---|
17 | * In this file, only the polar case is implemented. This corresponds to
|
---|
18 | * EPSG:9810 (Polar Stereographic Variant A) and EPSG:9829 (Polar Stereographic
|
---|
19 | * Variant B).
|
---|
20 | * <p>
|
---|
21 | *
|
---|
22 | * It is required, that the latitude of natural origin has the value +/-90 degrees.
|
---|
23 | * <p>
|
---|
24 | *
|
---|
25 | * This class has been derived from the implementation of the Geotools project;
|
---|
26 | * git 8cbf52d, org.geotools.referencing.operation.projection.PolarStereographic
|
---|
27 | * at the time of migration.
|
---|
28 | * <p>
|
---|
29 | *
|
---|
30 | * <b>References:</b>
|
---|
31 | * <ul>
|
---|
32 | * <li>John P. Snyder (Map Projections - A Working Manual,<br>
|
---|
33 | * U.S. Geological Survey Professional Paper 1395, 1987)</li>
|
---|
34 | * <li>"Coordinate Conversions and Transformations including Formulas",<br>
|
---|
35 | * EPSG Guidence Note Number 7, Version 19.</li>
|
---|
36 | * <li>Gerald Evenden. <A HREF="http://members.bellatlantic.net/~vze2hc4d/proj4/sterea.pdf">
|
---|
37 | * "Supplementary PROJ.4 Notes - Oblique Stereographic Alternative"</A></li>
|
---|
38 | * <li>Krakiwsky, E.J., D.B. Thomson, and R.R. Steeves. 1977. A Manual
|
---|
39 | * For Geodetic Coordinate Transformations in the Maritimes.
|
---|
40 | * Geodesy and Geomatics Engineering, UNB. Technical Report No. 48.</li>
|
---|
41 | * <li>Thomson, D.B., M.P. Mepham and R.R. Steeves. 1977.
|
---|
42 | * The Stereographic Double Projection.
|
---|
43 | * Geodesy and Geomatics Engineereng, UNB. Technical Report No. 46.</li>
|
---|
44 | * </ul>
|
---|
45 | *
|
---|
46 | * @author André Gosselin
|
---|
47 | * @author Martin Desruisseaux (PMO, IRD)
|
---|
48 | * @author Rueben Schulz
|
---|
49 | *
|
---|
50 | * @see <A HREF="http://mathworld.wolfram.com/StereographicProjection.html">Stereographic projection on MathWorld</A>
|
---|
51 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/polar_stereographic.html">Polar_Stereographic</A>
|
---|
52 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/oblique_stereographic.html">Oblique_Stereographic</A>
|
---|
53 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/stereographic.html">Stereographic</A>
|
---|
54 | * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/random_issues.html#stereographic">Some Random Stereographic Issues</A>
|
---|
55 | *
|
---|
56 | * @see DoubleStereographic
|
---|
57 | * @since 9419
|
---|
58 | */
|
---|
59 | public class PolarStereographic extends AbstractProj {
|
---|
60 | /**
|
---|
61 | * Maximum number of iterations for iterative computations.
|
---|
62 | */
|
---|
63 | private static final int MAXIMUM_ITERATIONS = 15;
|
---|
64 |
|
---|
65 | /**
|
---|
66 | * Difference allowed in iterative computations.
|
---|
67 | */
|
---|
68 | private static final double ITERATION_TOLERANCE = 1E-10;
|
---|
69 |
|
---|
70 | /**
|
---|
71 | * Maximum difference allowed when comparing real numbers.
|
---|
72 | */
|
---|
73 | private static final double EPSILON = 1E-8;
|
---|
74 |
|
---|
75 | /**
|
---|
76 | * A constant used in the transformations.
|
---|
77 | */
|
---|
78 | private double k0;
|
---|
79 |
|
---|
80 | /**
|
---|
81 | * {@code true} if this projection is for the south pole, or {@code false}
|
---|
82 | * if it is for the north pole.
|
---|
83 | */
|
---|
84 | boolean southPole;
|
---|
85 |
|
---|
86 | @Override
|
---|
87 | public String getName() {
|
---|
88 | return tr("Polar Stereographic");
|
---|
89 | }
|
---|
90 |
|
---|
91 | @Override
|
---|
92 | public String getProj4Id() {
|
---|
93 | return "stere";
|
---|
94 | }
|
---|
95 |
|
---|
96 | @Override
|
---|
97 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
98 | super.initialize(params);
|
---|
99 | if (params.lat0 == null)
|
---|
100 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
|
---|
101 | if (params.lat0 != 90.0 && params.lat0 != -90.0)
|
---|
102 | throw new ProjectionConfigurationException(
|
---|
103 | tr("Polar Stereographic: Parameter ''{0}'' must be 90 or -90.", "lat_0"));
|
---|
104 | // Latitude of true scale, in radians;
|
---|
105 | double latitudeTrueScale;
|
---|
106 | if (params.lat_ts == null) {
|
---|
107 | latitudeTrueScale = (params.lat0 < 0) ? -Math.PI/2 : Math.PI/2;
|
---|
108 | } else {
|
---|
109 | latitudeTrueScale = Math.toRadians(params.lat_ts);
|
---|
110 | }
|
---|
111 | southPole = latitudeTrueScale < 0;
|
---|
112 |
|
---|
113 | // Computes coefficients.
|
---|
114 | double latitudeTrueScaleAbs = Math.abs(latitudeTrueScale);
|
---|
115 | if (Math.abs(latitudeTrueScaleAbs - Math.PI/2) >= EPSILON) {
|
---|
116 | final double t = Math.sin(latitudeTrueScaleAbs);
|
---|
117 | k0 = msfn(t, Math.cos(latitudeTrueScaleAbs)) /
|
---|
118 | tsfn(latitudeTrueScaleAbs, t); // Derives from (21-32 and 21-33)
|
---|
119 | } else {
|
---|
120 | // True scale at pole (part of (21-33))
|
---|
121 | k0 = 2.0 / Math.sqrt(Math.pow(1+e, 1+e)*
|
---|
122 | Math.pow(1-e, 1-e));
|
---|
123 | }
|
---|
124 | }
|
---|
125 |
|
---|
126 | @Override
|
---|
127 | public double[] project(double y, double x) {
|
---|
128 | final double sinlat = Math.sin(y);
|
---|
129 | final double coslon = Math.cos(x);
|
---|
130 | final double sinlon = Math.sin(x);
|
---|
131 | if (southPole) {
|
---|
132 | final double rho = k0 * tsfn(-y, -sinlat);
|
---|
133 | x = rho * sinlon;
|
---|
134 | y = rho * coslon;
|
---|
135 | } else {
|
---|
136 | final double rho = k0 * tsfn(y, sinlat);
|
---|
137 | x = rho * sinlon;
|
---|
138 | y = -rho * coslon;
|
---|
139 | }
|
---|
140 | return new double[] {x, y};
|
---|
141 | }
|
---|
142 |
|
---|
143 | @Override
|
---|
144 | public double[] invproject(double x, double y) {
|
---|
145 | final double rho = Math.hypot(x, y);
|
---|
146 | if (southPole) {
|
---|
147 | y = -y;
|
---|
148 | }
|
---|
149 | /*
|
---|
150 | * Compute latitude using iterative technique.
|
---|
151 | */
|
---|
152 | final double t = rho/k0;
|
---|
153 | final double halfe = e/2.0;
|
---|
154 | double phi0 = 0;
|
---|
155 | for (int i = MAXIMUM_ITERATIONS;;) {
|
---|
156 | final double esinphi = e * Math.sin(phi0);
|
---|
157 | final double phi = (Math.PI/2) - 2.0*Math.atan(t*Math.pow((1-esinphi)/(1+esinphi), halfe));
|
---|
158 | if (Math.abs(phi-phi0) < ITERATION_TOLERANCE) {
|
---|
159 | x = (Math.abs(rho) < EPSILON) ? 0.0 : Math.atan2(x, -y);
|
---|
160 | y = southPole ? -phi : phi;
|
---|
161 | break;
|
---|
162 | }
|
---|
163 | phi0 = phi;
|
---|
164 | if (--i < 0) {
|
---|
165 | throw new IllegalStateException("no convergence for x="+x+", y="+y);
|
---|
166 | }
|
---|
167 | }
|
---|
168 | return new double[] {y, x};
|
---|
169 | }
|
---|
170 |
|
---|
171 | @Override
|
---|
172 | public Bounds getAlgorithmBounds() {
|
---|
173 | final double cut = 60;
|
---|
174 | if (southPole) {
|
---|
175 | return new Bounds(-90, -180, cut, 180, false);
|
---|
176 | } else {
|
---|
177 | return new Bounds(-cut, -180, 90, 180, false);
|
---|
178 | }
|
---|
179 | }
|
---|
180 | }
|
---|