source: josm/trunk/src/org/openstreetmap/josm/data/projection/UTM.java@ 2103

Last change on this file since 2103 was 2017, checked in by Gubaer, 15 years ago

removed OptionPaneUtil
cleanup of deprecated Layer API
cleanup of deprecated APIs in OsmPrimitive and Way
cleanup of imports

File size: 12.7 KB
Line 
1// License: GPL. Copyright 2007 by Immanuel Scholz and others
2package org.openstreetmap.josm.data.projection;
3
4import static org.openstreetmap.josm.tools.I18n.tr;
5
6import org.openstreetmap.josm.data.Bounds;
7import org.openstreetmap.josm.data.coor.EastNorth;
8import org.openstreetmap.josm.data.coor.LatLon;
9
10/**
11 * Directly use latitude / longitude values as x/y.
12 *
13 * @author Dirk Stöcker
14 * code based on JavaScript from Chuck Taylor
15 */
16public class UTM implements Projection {
17
18 final private double UTMScaleFactor = 0.9996;
19
20 /* Ellipsoid model constants (WGS84) - TODO Use Elliposid class here too */
21 final private double sm_EccSquared = 6.69437999013e-03;
22
23 /*
24 * ArcLengthOfMeridian
25 *
26 * Computes the ellipsoidal distance from the equator to a point at a
27 * given latitude.
28 *
29 * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
30 * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
31 *
32 * Inputs:
33 * phi - Latitude of the point, in radians.
34 *
35 * Globals:
36 * Ellipsoid.GRS80.a - Ellipsoid model major axis.
37 * Ellipsoid.GRS80.b - Ellipsoid model minor axis.
38 *
39 * Returns:
40 * The ellipsoidal distance of the point from the equator, in meters.
41 *
42 */
43 private double ArcLengthOfMeridian(double phi)
44 {
45 /* Precalculate n */
46 double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
47
48 /* Precalculate alpha */
49 double alpha = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
50 * (1.0 + (Math.pow (n, 2.0) / 4.0) + (Math.pow (n, 4.0) / 64.0));
51
52 /* Precalculate beta */
53 double beta = (-3.0 * n / 2.0) + (9.0 * Math.pow (n, 3.0) / 16.0)
54 + (-3.0 * Math.pow (n, 5.0) / 32.0);
55
56 /* Precalculate gamma */
57 double gamma = (15.0 * Math.pow (n, 2.0) / 16.0)
58 + (-15.0 * Math.pow (n, 4.0) / 32.0);
59
60 /* Precalculate delta */
61 double delta = (-35.0 * Math.pow (n, 3.0) / 48.0)
62 + (105.0 * Math.pow (n, 5.0) / 256.0);
63
64 /* Precalculate epsilon */
65 double epsilon = (315.0 * Math.pow (n, 4.0) / 512.0);
66
67 /* Now calculate the sum of the series and return */
68 return alpha
69 * (phi + (beta * Math.sin (2.0 * phi))
70 + (gamma * Math.sin (4.0 * phi))
71 + (delta * Math.sin (6.0 * phi))
72 + (epsilon * Math.sin (8.0 * phi)));
73 }
74
75 /*
76 * UTMCentralMeridian
77 *
78 * Determines the central meridian for the given UTM zone.
79 *
80 * Inputs:
81 * zone - An integer value designating the UTM zone, range [1,60].
82 *
83 * Returns:
84 * The central meridian for the given UTM zone, in radians, or zero
85 * if the UTM zone parameter is outside the range [1,60].
86 * Range of the central meridian is the radian equivalent of [-177,+177].
87 *
88 */
89 private double UTMCentralMeridian(int zone)
90 {
91 return Math.toRadians(-183.0 + (zone * 6.0));
92 }
93 private double UTMCentralMeridianDeg(int zone)
94 {
95 return -183.0 + (zone * 6.0);
96 }
97
98 /*
99 * FootpointLatitude
100 *
101 * Computes the footpoint latitude for use in converting transverse
102 * Mercator coordinates to ellipsoidal coordinates.
103 *
104 * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
105 * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
106 *
107 * Inputs:
108 * y - The UTM northing coordinate, in meters.
109 *
110 * Returns:
111 * The footpoint latitude, in radians.
112 *
113 */
114 private double FootpointLatitude(double y)
115 {
116 /* Precalculate n (Eq. 10.18) */
117 double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
118
119 /* Precalculate alpha_ (Eq. 10.22) */
120 /* (Same as alpha in Eq. 10.17) */
121 double alpha_ = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
122 * (1 + (Math.pow (n, 2.0) / 4) + (Math.pow (n, 4.0) / 64));
123
124 /* Precalculate y_ (Eq. 10.23) */
125 double y_ = y / alpha_;
126
127 /* Precalculate beta_ (Eq. 10.22) */
128 double beta_ = (3.0 * n / 2.0) + (-27.0 * Math.pow (n, 3.0) / 32.0)
129 + (269.0 * Math.pow (n, 5.0) / 512.0);
130
131 /* Precalculate gamma_ (Eq. 10.22) */
132 double gamma_ = (21.0 * Math.pow (n, 2.0) / 16.0)
133 + (-55.0 * Math.pow (n, 4.0) / 32.0);
134
135 /* Precalculate delta_ (Eq. 10.22) */
136 double delta_ = (151.0 * Math.pow (n, 3.0) / 96.0)
137 + (-417.0 * Math.pow (n, 5.0) / 128.0);
138
139 /* Precalculate epsilon_ (Eq. 10.22) */
140 double epsilon_ = (1097.0 * Math.pow (n, 4.0) / 512.0);
141
142 /* Now calculate the sum of the series (Eq. 10.21) */
143 return y_ + (beta_ * Math.sin (2.0 * y_))
144 + (gamma_ * Math.sin (4.0 * y_))
145 + (delta_ * Math.sin (6.0 * y_))
146 + (epsilon_ * Math.sin (8.0 * y_));
147 }
148
149 /*
150 * MapLatLonToXY
151 *
152 * Converts a latitude/longitude pair to x and y coordinates in the
153 * Transverse Mercator projection. Note that Transverse Mercator is not
154 * the same as UTM; a scale factor is required to convert between them.
155 *
156 * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
157 * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
158 *
159 * Inputs:
160 * phi - Latitude of the point, in radians.
161 * lambda - Longitude of the point, in radians.
162 * lambda0 - Longitude of the central meridian to be used, in radians.
163 *
164 * Outputs:
165 * xy - A 2-element array containing the x and y coordinates
166 * of the computed point.
167 *
168 * Returns:
169 * The function does not return a value.
170 *
171 */
172 public EastNorth MapLatLonToXY(double phi, double lambda, double lambda0)
173 {
174 /* Precalculate ep2 */
175 double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0)) / Math.pow (Ellipsoid.GRS80.b, 2.0);
176
177 /* Precalculate nu2 */
178 double nu2 = ep2 * Math.pow (Math.cos (phi), 2.0);
179
180 /* Precalculate N */
181 double N = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nu2));
182
183 /* Precalculate t */
184 double t = Math.tan (phi);
185 double t2 = t * t;
186 double tmp = (t2 * t2 * t2) - Math.pow (t, 6.0);
187
188 /* Precalculate l */
189 double l = lambda - lambda0;
190
191 /* Precalculate coefficients for l**n in the equations below
192 so a normal human being can read the expressions for easting
193 and northing
194 -- l**1 and l**2 have coefficients of 1.0 */
195 double l3coef = 1.0 - t2 + nu2;
196
197 double l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2);
198
199 double l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2
200 - 58.0 * t2 * nu2;
201
202 double l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2
203 - 330.0 * t2 * nu2;
204
205 double l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2);
206
207 double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
208
209 return new EastNorth(
210 /* Calculate easting (x) */
211 N * Math.cos (phi) * l
212 + (N / 6.0 * Math.pow (Math.cos (phi), 3.0) * l3coef * Math.pow (l, 3.0))
213 + (N / 120.0 * Math.pow (Math.cos (phi), 5.0) * l5coef * Math.pow (l, 5.0))
214 + (N / 5040.0 * Math.pow (Math.cos (phi), 7.0) * l7coef * Math.pow (l, 7.0)),
215 /* Calculate northing (y) */
216 ArcLengthOfMeridian (phi)
217 + (t / 2.0 * N * Math.pow (Math.cos (phi), 2.0) * Math.pow (l, 2.0))
218 + (t / 24.0 * N * Math.pow (Math.cos (phi), 4.0) * l4coef * Math.pow (l, 4.0))
219 + (t / 720.0 * N * Math.pow (Math.cos (phi), 6.0) * l6coef * Math.pow (l, 6.0))
220 + (t / 40320.0 * N * Math.pow (Math.cos (phi), 8.0) * l8coef * Math.pow (l, 8.0)));
221 }
222
223 /*
224 * MapXYToLatLon
225 *
226 * Converts x and y coordinates in the Transverse Mercator projection to
227 * a latitude/longitude pair. Note that Transverse Mercator is not
228 * the same as UTM; a scale factor is required to convert between them.
229 *
230 * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
231 * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
232 *
233 * Inputs:
234 * x - The easting of the point, in meters.
235 * y - The northing of the point, in meters.
236 * lambda0 - Longitude of the central meridian to be used, in radians.
237 *
238 * Outputs:
239 * philambda - A 2-element containing the latitude and longitude
240 * in radians.
241 *
242 * Returns:
243 * The function does not return a value.
244 *
245 * Remarks:
246 * The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
247 * N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect
248 * to the footpoint latitude phif.
249 *
250 * x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
251 * to optimize computations.
252 *
253 */
254 public LatLon MapXYToLatLon(double x, double y, double lambda0)
255 {
256 /* Get the value of phif, the footpoint latitude. */
257 double phif = FootpointLatitude (y);
258
259 /* Precalculate ep2 */
260 double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0))
261 / Math.pow (Ellipsoid.GRS80.b, 2.0);
262
263 /* Precalculate cos (phif) */
264 double cf = Math.cos (phif);
265
266 /* Precalculate nuf2 */
267 double nuf2 = ep2 * Math.pow (cf, 2.0);
268
269 /* Precalculate Nf and initialize Nfpow */
270 double Nf = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nuf2));
271 double Nfpow = Nf;
272
273 /* Precalculate tf */
274 double tf = Math.tan (phif);
275 double tf2 = tf * tf;
276 double tf4 = tf2 * tf2;
277
278 /* Precalculate fractional coefficients for x**n in the equations
279 below to simplify the expressions for latitude and longitude. */
280 double x1frac = 1.0 / (Nfpow * cf);
281
282 Nfpow *= Nf; /* now equals Nf**2) */
283 double x2frac = tf / (2.0 * Nfpow);
284
285 Nfpow *= Nf; /* now equals Nf**3) */
286 double x3frac = 1.0 / (6.0 * Nfpow * cf);
287
288 Nfpow *= Nf; /* now equals Nf**4) */
289 double x4frac = tf / (24.0 * Nfpow);
290
291 Nfpow *= Nf; /* now equals Nf**5) */
292 double x5frac = 1.0 / (120.0 * Nfpow * cf);
293
294 Nfpow *= Nf; /* now equals Nf**6) */
295 double x6frac = tf / (720.0 * Nfpow);
296
297 Nfpow *= Nf; /* now equals Nf**7) */
298 double x7frac = 1.0 / (5040.0 * Nfpow * cf);
299
300 Nfpow *= Nf; /* now equals Nf**8) */
301 double x8frac = tf / (40320.0 * Nfpow);
302
303 /* Precalculate polynomial coefficients for x**n.
304 -- x**1 does not have a polynomial coefficient. */
305 double x2poly = -1.0 - nuf2;
306 double x3poly = -1.0 - 2 * tf2 - nuf2;
307 double x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2);
308 double x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2;
309 double x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2;
310 double x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2);
311 double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2);
312
313 return new LatLon(
314 /* Calculate latitude */
315 Math.toDegrees(
316 phif + x2frac * x2poly * (x * x)
317 + x4frac * x4poly * Math.pow (x, 4.0)
318 + x6frac * x6poly * Math.pow (x, 6.0)
319 + x8frac * x8poly * Math.pow (x, 8.0)),
320 Math.toDegrees(
321 /* Calculate longitude */
322 lambda0 + x1frac * x
323 + x3frac * x3poly * Math.pow (x, 3.0)
324 + x5frac * x5poly * Math.pow (x, 5.0)
325 + x7frac * x7poly * Math.pow (x, 7.0)));
326 }
327
328 public EastNorth latlon2eastNorth(LatLon p) {
329 EastNorth a = MapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), UTMCentralMeridian(getzone()));
330 return new EastNorth(a.east() * UTMScaleFactor + 3500000.0, a.north() * UTMScaleFactor);
331 }
332
333 public LatLon eastNorth2latlon(EastNorth p) {
334 return MapXYToLatLon((p.east()-3500000.0)/UTMScaleFactor, p.north()/UTMScaleFactor, UTMCentralMeridian(getzone()));
335 }
336
337 @Override public String toString() {
338 return tr("UTM Zone {0}", getzone());
339 }
340
341 /* TODO - support all UTM's not only zone 33 */
342 public int getzone()
343 {
344 return 33;
345 }
346
347 public String toCode() {
348 return "EPSG:325833";
349 }
350
351 public String getCacheDirectoryName() {
352 return "epsg325833";
353 }
354
355 public Bounds getWorldBoundsLatLon()
356 {
357 return new Bounds(
358 new LatLon(-85.0, UTMCentralMeridianDeg(getzone())-5.0),
359 new LatLon(85.0, UTMCentralMeridianDeg(getzone())+5.0));
360 }
361}
Note: See TracBrowser for help on using the repository browser.