1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import static java.lang.Math.cos;
|
---|
5 | import static java.lang.Math.pow;
|
---|
6 | import static java.lang.Math.sin;
|
---|
7 | import static java.lang.Math.sqrt;
|
---|
8 | import static java.lang.Math.tan;
|
---|
9 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
10 |
|
---|
11 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
12 |
|
---|
13 | /**
|
---|
14 | * Transverse Mercator projection.
|
---|
15 | *
|
---|
16 | * @author Dirk Stöcker
|
---|
17 | * code based on JavaScript from Chuck Taylor
|
---|
18 | *
|
---|
19 | */
|
---|
20 | public class TransverseMercator implements Proj {
|
---|
21 |
|
---|
22 | protected double a, b;
|
---|
23 |
|
---|
24 | @Override
|
---|
25 | public String getName() {
|
---|
26 | return tr("Transverse Mercator");
|
---|
27 | }
|
---|
28 |
|
---|
29 | @Override
|
---|
30 | public String getProj4Id() {
|
---|
31 | return "tmerc";
|
---|
32 | }
|
---|
33 |
|
---|
34 | @Override
|
---|
35 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
36 | this.a = params.ellps.a;
|
---|
37 | this.b = params.ellps.b;
|
---|
38 | }
|
---|
39 |
|
---|
40 | /**
|
---|
41 | * Converts a latitude/longitude pair to x and y coordinates in the
|
---|
42 | * Transverse Mercator projection. Note that Transverse Mercator is not
|
---|
43 | * the same as UTM; a scale factor is required to convert between them.
|
---|
44 | *
|
---|
45 | * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
|
---|
46 | * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
|
---|
47 | *
|
---|
48 | * @param phi Latitude of the point, in radians
|
---|
49 | * @param lambda Longitude of the point, in radians
|
---|
50 | * @return A 2-element array containing the x and y coordinates
|
---|
51 | * of the computed point
|
---|
52 | */
|
---|
53 | @Override
|
---|
54 | public double[] project(double phi, double lambda) {
|
---|
55 |
|
---|
56 | /* Precalculate ep2 */
|
---|
57 | double ep2 = (pow(a, 2.0) - pow(b, 2.0)) / pow(b, 2.0);
|
---|
58 |
|
---|
59 | /* Precalculate nu2 */
|
---|
60 | double nu2 = ep2 * pow(cos(phi), 2.0);
|
---|
61 |
|
---|
62 | /* Precalculate N / a */
|
---|
63 | double N_a = a / (b * sqrt(1 + nu2));
|
---|
64 |
|
---|
65 | /* Precalculate t */
|
---|
66 | double t = tan(phi);
|
---|
67 | double t2 = t * t;
|
---|
68 |
|
---|
69 | /* Precalculate l */
|
---|
70 | double l = lambda;
|
---|
71 |
|
---|
72 | /* Precalculate coefficients for l**n in the equations below
|
---|
73 | so a normal human being can read the expressions for easting
|
---|
74 | and northing
|
---|
75 | -- l**1 and l**2 have coefficients of 1.0 */
|
---|
76 | double l3coef = 1.0 - t2 + nu2;
|
---|
77 |
|
---|
78 | double l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2);
|
---|
79 |
|
---|
80 | double l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2
|
---|
81 | - 58.0 * t2 * nu2;
|
---|
82 |
|
---|
83 | double l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2
|
---|
84 | - 330.0 * t2 * nu2;
|
---|
85 |
|
---|
86 | double l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2);
|
---|
87 |
|
---|
88 | double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
|
---|
89 |
|
---|
90 | return new double[] {
|
---|
91 | /* Calculate easting (x) */
|
---|
92 | N_a * cos(phi) * l
|
---|
93 | + (N_a / 6.0 * pow(cos(phi), 3.0) * l3coef * pow(l, 3.0))
|
---|
94 | + (N_a / 120.0 * pow(cos(phi), 5.0) * l5coef * pow(l, 5.0))
|
---|
95 | + (N_a / 5040.0 * pow(cos(phi), 7.0) * l7coef * pow(l, 7.0)),
|
---|
96 | /* Calculate northing (y) */
|
---|
97 | ArcLengthOfMeridian (phi) / a
|
---|
98 | + (t / 2.0 * N_a * pow(cos(phi), 2.0) * pow(l, 2.0))
|
---|
99 | + (t / 24.0 * N_a * pow(cos(phi), 4.0) * l4coef * pow(l, 4.0))
|
---|
100 | + (t / 720.0 * N_a * pow(cos(phi), 6.0) * l6coef * pow(l, 6.0))
|
---|
101 | + (t / 40320.0 * N_a * pow(cos(phi), 8.0) * l8coef * pow(l, 8.0)) };
|
---|
102 | }
|
---|
103 |
|
---|
104 | /**
|
---|
105 | * Converts x and y coordinates in the Transverse Mercator projection to
|
---|
106 | * a latitude/longitude pair. Note that Transverse Mercator is not
|
---|
107 | * the same as UTM; a scale factor is required to convert between them.
|
---|
108 | *
|
---|
109 | * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
|
---|
110 | * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
|
---|
111 | *
|
---|
112 | * Remarks:
|
---|
113 | * The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
|
---|
114 | * N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect
|
---|
115 | * to the footpoint latitude phif.
|
---|
116 | *
|
---|
117 | * x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
|
---|
118 | * to optimize computations.
|
---|
119 | *
|
---|
120 | * @param x The easting of the point, in meters, divided by the semi major axis of the ellipsoid
|
---|
121 | * @param y The northing of the point, in meters, divided by the semi major axis of the ellipsoid
|
---|
122 | * @return A 2-element containing the latitude and longitude
|
---|
123 | * in radians
|
---|
124 | */
|
---|
125 | @Override
|
---|
126 | public double[] invproject(double x, double y) {
|
---|
127 | /* Get the value of phif, the footpoint latitude. */
|
---|
128 | double phif = footpointLatitude(y);
|
---|
129 |
|
---|
130 | /* Precalculate ep2 */
|
---|
131 | double ep2 = (a*a - b*b)
|
---|
132 | / (b*b);
|
---|
133 |
|
---|
134 | /* Precalculate cos(phif) */
|
---|
135 | double cf = cos(phif);
|
---|
136 |
|
---|
137 | /* Precalculate nuf2 */
|
---|
138 | double nuf2 = ep2 * pow(cf, 2.0);
|
---|
139 |
|
---|
140 | /* Precalculate Nf / a and initialize Nfpow */
|
---|
141 | double Nf_a = a / (b * sqrt(1 + nuf2));
|
---|
142 | double Nfpow = Nf_a;
|
---|
143 |
|
---|
144 | /* Precalculate tf */
|
---|
145 | double tf = tan(phif);
|
---|
146 | double tf2 = tf * tf;
|
---|
147 | double tf4 = tf2 * tf2;
|
---|
148 |
|
---|
149 | /* Precalculate fractional coefficients for x**n in the equations
|
---|
150 | below to simplify the expressions for latitude and longitude. */
|
---|
151 | double x1frac = 1.0 / (Nfpow * cf);
|
---|
152 |
|
---|
153 | Nfpow *= Nf_a; /* now equals Nf**2) */
|
---|
154 | double x2frac = tf / (2.0 * Nfpow);
|
---|
155 |
|
---|
156 | Nfpow *= Nf_a; /* now equals Nf**3) */
|
---|
157 | double x3frac = 1.0 / (6.0 * Nfpow * cf);
|
---|
158 |
|
---|
159 | Nfpow *= Nf_a; /* now equals Nf**4) */
|
---|
160 | double x4frac = tf / (24.0 * Nfpow);
|
---|
161 |
|
---|
162 | Nfpow *= Nf_a; /* now equals Nf**5) */
|
---|
163 | double x5frac = 1.0 / (120.0 * Nfpow * cf);
|
---|
164 |
|
---|
165 | Nfpow *= Nf_a; /* now equals Nf**6) */
|
---|
166 | double x6frac = tf / (720.0 * Nfpow);
|
---|
167 |
|
---|
168 | Nfpow *= Nf_a; /* now equals Nf**7) */
|
---|
169 | double x7frac = 1.0 / (5040.0 * Nfpow * cf);
|
---|
170 |
|
---|
171 | Nfpow *= Nf_a; /* now equals Nf**8) */
|
---|
172 | double x8frac = tf / (40320.0 * Nfpow);
|
---|
173 |
|
---|
174 | /* Precalculate polynomial coefficients for x**n.
|
---|
175 | -- x**1 does not have a polynomial coefficient. */
|
---|
176 | double x2poly = -1.0 - nuf2;
|
---|
177 | double x3poly = -1.0 - 2 * tf2 - nuf2;
|
---|
178 | double x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2);
|
---|
179 | double x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2;
|
---|
180 | double x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2;
|
---|
181 | double x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2);
|
---|
182 | double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2);
|
---|
183 |
|
---|
184 | return new double[] {
|
---|
185 | /* Calculate latitude */
|
---|
186 | phif + x2frac * x2poly * (x * x)
|
---|
187 | + x4frac * x4poly * pow(x, 4.0)
|
---|
188 | + x6frac * x6poly * pow(x, 6.0)
|
---|
189 | + x8frac * x8poly * pow(x, 8.0),
|
---|
190 | /* Calculate longitude */
|
---|
191 | x1frac * x
|
---|
192 | + x3frac * x3poly * pow(x, 3.0)
|
---|
193 | + x5frac * x5poly * pow(x, 5.0)
|
---|
194 | + x7frac * x7poly * pow(x, 7.0) };
|
---|
195 | }
|
---|
196 |
|
---|
197 | /**
|
---|
198 | * ArcLengthOfMeridian
|
---|
199 | *
|
---|
200 | * Computes the ellipsoidal distance from the equator to a point at a
|
---|
201 | * given latitude.
|
---|
202 | *
|
---|
203 | * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
|
---|
204 | * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
|
---|
205 | *
|
---|
206 | * @param phi Latitude of the point, in radians
|
---|
207 | * @return The ellipsoidal distance of the point from the equator
|
---|
208 | * (in meters, divided by the semi major axis of the ellipsoid)
|
---|
209 | */
|
---|
210 | private double ArcLengthOfMeridian(double phi) {
|
---|
211 | /* Precalculate n */
|
---|
212 | double n = (a - b) / (a + b);
|
---|
213 |
|
---|
214 | /* Precalculate alpha */
|
---|
215 | double alpha = ((a + b) / 2.0)
|
---|
216 | * (1.0 + (pow(n, 2.0) / 4.0) + (pow(n, 4.0) / 64.0));
|
---|
217 |
|
---|
218 | /* Precalculate beta */
|
---|
219 | double beta = (-3.0 * n / 2.0) + (9.0 * pow(n, 3.0) / 16.0)
|
---|
220 | + (-3.0 * pow(n, 5.0) / 32.0);
|
---|
221 |
|
---|
222 | /* Precalculate gamma */
|
---|
223 | double gamma = (15.0 * pow(n, 2.0) / 16.0)
|
---|
224 | + (-15.0 * pow(n, 4.0) / 32.0);
|
---|
225 |
|
---|
226 | /* Precalculate delta */
|
---|
227 | double delta = (-35.0 * pow(n, 3.0) / 48.0)
|
---|
228 | + (105.0 * pow(n, 5.0) / 256.0);
|
---|
229 |
|
---|
230 | /* Precalculate epsilon */
|
---|
231 | double epsilon = 315.0 * pow(n, 4.0) / 512.0;
|
---|
232 |
|
---|
233 | /* Now calculate the sum of the series and return */
|
---|
234 | return alpha
|
---|
235 | * (phi + (beta * sin(2.0 * phi))
|
---|
236 | + (gamma * sin(4.0 * phi))
|
---|
237 | + (delta * sin(6.0 * phi))
|
---|
238 | + (epsilon * sin(8.0 * phi)));
|
---|
239 | }
|
---|
240 |
|
---|
241 | /**
|
---|
242 | * FootpointLatitude
|
---|
243 | *
|
---|
244 | * Computes the footpoint latitude for use in converting transverse
|
---|
245 | * Mercator coordinates to ellipsoidal coordinates.
|
---|
246 | *
|
---|
247 | * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
|
---|
248 | * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
|
---|
249 | *
|
---|
250 | * @param y northing coordinate, in meters, divided by the semi major axis of the ellipsoid
|
---|
251 | * @return The footpoint latitude, in radians
|
---|
252 | */
|
---|
253 | private double footpointLatitude(double y) {
|
---|
254 | /* Precalculate n (Eq. 10.18) */
|
---|
255 | double n = (a - b) / (a + b);
|
---|
256 |
|
---|
257 | /* Precalculate alpha_ (Eq. 10.22) */
|
---|
258 | /* (Same as alpha in Eq. 10.17) */
|
---|
259 | double alpha_ = ((a + b) / 2.0)
|
---|
260 | * (1 + (pow(n, 2.0) / 4) + (pow(n, 4.0) / 64));
|
---|
261 |
|
---|
262 | /* Precalculate y_ (Eq. 10.23) */
|
---|
263 | double y_ = y / alpha_ * a;
|
---|
264 |
|
---|
265 | /* Precalculate beta_ (Eq. 10.22) */
|
---|
266 | double beta_ = (3.0 * n / 2.0) + (-27.0 * pow(n, 3.0) / 32.0)
|
---|
267 | + (269.0 * pow(n, 5.0) / 512.0);
|
---|
268 |
|
---|
269 | /* Precalculate gamma_ (Eq. 10.22) */
|
---|
270 | double gamma_ = (21.0 * pow(n, 2.0) / 16.0)
|
---|
271 | + (-55.0 * pow(n, 4.0) / 32.0);
|
---|
272 |
|
---|
273 | /* Precalculate delta_ (Eq. 10.22) */
|
---|
274 | double delta_ = (151.0 * pow(n, 3.0) / 96.0)
|
---|
275 | + (-417.0 * pow(n, 5.0) / 128.0);
|
---|
276 |
|
---|
277 | /* Precalculate epsilon_ (Eq. 10.22) */
|
---|
278 | double epsilon_ = 1097.0 * pow(n, 4.0) / 512.0;
|
---|
279 |
|
---|
280 | /* Now calculate the sum of the series (Eq. 10.21) */
|
---|
281 | return y_ + (beta_ * sin(2.0 * y_))
|
---|
282 | + (gamma_ * sin(4.0 * y_))
|
---|
283 | + (delta_ * sin(6.0 * y_))
|
---|
284 | + (epsilon_ * sin(8.0 * y_));
|
---|
285 | }
|
---|
286 | }
|
---|