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

Last change on this file since 2516 was 2516, checked in by stoecker, 14 years ago

close #4015 - Zoomlevel changes whenever the preference dialog is closed
File size: 15.0 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 java.awt.GridBagLayout;
7import java.util.Collection;
8import java.util.Collections;
9
10import javax.swing.JComboBox;
11import javax.swing.JLabel;
12import javax.swing.JPanel;
13
14import org.openstreetmap.josm.Main;
15import org.openstreetmap.josm.data.Bounds;
16import org.openstreetmap.josm.data.coor.EastNorth;
17import org.openstreetmap.josm.data.coor.LatLon;
18import org.openstreetmap.josm.tools.GBC;
19
20/**
21 * Directly use latitude / longitude values as x/y.
22 *
23 * @author Dirk Stöcker
24 * code based on JavaScript from Chuck Taylor
25 */
26public class UTM implements Projection, ProjectionSubPrefs {
27
28 public static final int DEFAULT_ZONE = 30;
29 private int zone = DEFAULT_ZONE;
30
31 final private double UTMScaleFactor = 0.9996;
32
33 /* Ellipsoid model constants (WGS84) - TODO Use Elliposid class here too */
34 final private double sm_EccSquared = 6.69437999013e-03;
35
36 /*
37 * ArcLengthOfMeridian
38 *
39 * Computes the ellipsoidal distance from the equator to a point at a
40 * given latitude.
41 *
42 * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
43 * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
44 *
45 * Inputs:
46 * phi - Latitude of the point, in radians.
47 *
48 * Globals:
49 * Ellipsoid.GRS80.a - Ellipsoid model major axis.
50 * Ellipsoid.GRS80.b - Ellipsoid model minor axis.
51 *
52 * Returns:
53 * The ellipsoidal distance of the point from the equator, in meters.
54 *
55 */
56 private double ArcLengthOfMeridian(double phi)
57 {
58 /* Precalculate n */
59 double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
60
61 /* Precalculate alpha */
62 double alpha = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
63 * (1.0 + (Math.pow (n, 2.0) / 4.0) + (Math.pow (n, 4.0) / 64.0));
64
65 /* Precalculate beta */
66 double beta = (-3.0 * n / 2.0) + (9.0 * Math.pow (n, 3.0) / 16.0)
67 + (-3.0 * Math.pow (n, 5.0) / 32.0);
68
69 /* Precalculate gamma */
70 double gamma = (15.0 * Math.pow (n, 2.0) / 16.0)
71 + (-15.0 * Math.pow (n, 4.0) / 32.0);
72
73 /* Precalculate delta */
74 double delta = (-35.0 * Math.pow (n, 3.0) / 48.0)
75 + (105.0 * Math.pow (n, 5.0) / 256.0);
76
77 /* Precalculate epsilon */
78 double epsilon = (315.0 * Math.pow (n, 4.0) / 512.0);
79
80 /* Now calculate the sum of the series and return */
81 return alpha
82 * (phi + (beta * Math.sin (2.0 * phi))
83 + (gamma * Math.sin (4.0 * phi))
84 + (delta * Math.sin (6.0 * phi))
85 + (epsilon * Math.sin (8.0 * phi)));
86 }
87
88 /*
89 * UTMCentralMeridian
90 *
91 * Determines the central meridian for the given UTM zone.
92 *
93 * Inputs:
94 * zone - An integer value designating the UTM zone, range [1,60].
95 *
96 * Returns:
97 * The central meridian for the given UTM zone, in radians, or zero
98 * if the UTM zone parameter is outside the range [1,60].
99 * Range of the central meridian is the radian equivalent of [-177,+177].
100 *
101 */
102 private double UTMCentralMeridian(int zone)
103 {
104 return Math.toRadians(-183.0 + (zone * 6.0));
105 }
106 private double UTMCentralMeridianDeg(int zone)
107 {
108 return -183.0 + (zone * 6.0);
109 }
110
111 /*
112 * FootpointLatitude
113 *
114 * Computes the footpoint latitude for use in converting transverse
115 * Mercator coordinates to ellipsoidal coordinates.
116 *
117 * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
118 * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
119 *
120 * Inputs:
121 * y - The UTM northing coordinate, in meters.
122 *
123 * Returns:
124 * The footpoint latitude, in radians.
125 *
126 */
127 private double FootpointLatitude(double y)
128 {
129 /* Precalculate n (Eq. 10.18) */
130 double n = (Ellipsoid.GRS80.a - Ellipsoid.GRS80.b) / (Ellipsoid.GRS80.a + Ellipsoid.GRS80.b);
131
132 /* Precalculate alpha_ (Eq. 10.22) */
133 /* (Same as alpha in Eq. 10.17) */
134 double alpha_ = ((Ellipsoid.GRS80.a + Ellipsoid.GRS80.b) / 2.0)
135 * (1 + (Math.pow (n, 2.0) / 4) + (Math.pow (n, 4.0) / 64));
136
137 /* Precalculate y_ (Eq. 10.23) */
138 double y_ = y / alpha_;
139
140 /* Precalculate beta_ (Eq. 10.22) */
141 double beta_ = (3.0 * n / 2.0) + (-27.0 * Math.pow (n, 3.0) / 32.0)
142 + (269.0 * Math.pow (n, 5.0) / 512.0);
143
144 /* Precalculate gamma_ (Eq. 10.22) */
145 double gamma_ = (21.0 * Math.pow (n, 2.0) / 16.0)
146 + (-55.0 * Math.pow (n, 4.0) / 32.0);
147
148 /* Precalculate delta_ (Eq. 10.22) */
149 double delta_ = (151.0 * Math.pow (n, 3.0) / 96.0)
150 + (-417.0 * Math.pow (n, 5.0) / 128.0);
151
152 /* Precalculate epsilon_ (Eq. 10.22) */
153 double epsilon_ = (1097.0 * Math.pow (n, 4.0) / 512.0);
154
155 /* Now calculate the sum of the series (Eq. 10.21) */
156 return y_ + (beta_ * Math.sin (2.0 * y_))
157 + (gamma_ * Math.sin (4.0 * y_))
158 + (delta_ * Math.sin (6.0 * y_))
159 + (epsilon_ * Math.sin (8.0 * y_));
160 }
161
162 /*
163 * MapLatLonToXY
164 *
165 * Converts a latitude/longitude pair to x and y coordinates in the
166 * Transverse Mercator projection. Note that Transverse Mercator is not
167 * the same as UTM; a scale factor is required to convert between them.
168 *
169 * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
170 * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
171 *
172 * Inputs:
173 * phi - Latitude of the point, in radians.
174 * lambda - Longitude of the point, in radians.
175 * lambda0 - Longitude of the central meridian to be used, in radians.
176 *
177 * Outputs:
178 * xy - A 2-element array containing the x and y coordinates
179 * of the computed point.
180 *
181 * Returns:
182 * The function does not return a value.
183 *
184 */
185 public EastNorth MapLatLonToXY(double phi, double lambda, double lambda0)
186 {
187 /* Precalculate ep2 */
188 double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0)) / Math.pow (Ellipsoid.GRS80.b, 2.0);
189
190 /* Precalculate nu2 */
191 double nu2 = ep2 * Math.pow (Math.cos (phi), 2.0);
192
193 /* Precalculate N */
194 double N = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nu2));
195
196 /* Precalculate t */
197 double t = Math.tan (phi);
198 double t2 = t * t;
199 double tmp = (t2 * t2 * t2) - Math.pow (t, 6.0);
200
201 /* Precalculate l */
202 double l = lambda - lambda0;
203
204 /* Precalculate coefficients for l**n in the equations below
205 so a normal human being can read the expressions for easting
206 and northing
207 -- l**1 and l**2 have coefficients of 1.0 */
208 double l3coef = 1.0 - t2 + nu2;
209
210 double l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2);
211
212 double l5coef = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2
213 - 58.0 * t2 * nu2;
214
215 double l6coef = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2
216 - 330.0 * t2 * nu2;
217
218 double l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2);
219
220 double l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
221
222 return new EastNorth(
223 /* Calculate easting (x) */
224 N * Math.cos (phi) * l
225 + (N / 6.0 * Math.pow (Math.cos (phi), 3.0) * l3coef * Math.pow (l, 3.0))
226 + (N / 120.0 * Math.pow (Math.cos (phi), 5.0) * l5coef * Math.pow (l, 5.0))
227 + (N / 5040.0 * Math.pow (Math.cos (phi), 7.0) * l7coef * Math.pow (l, 7.0)),
228 /* Calculate northing (y) */
229 ArcLengthOfMeridian (phi)
230 + (t / 2.0 * N * Math.pow (Math.cos (phi), 2.0) * Math.pow (l, 2.0))
231 + (t / 24.0 * N * Math.pow (Math.cos (phi), 4.0) * l4coef * Math.pow (l, 4.0))
232 + (t / 720.0 * N * Math.pow (Math.cos (phi), 6.0) * l6coef * Math.pow (l, 6.0))
233 + (t / 40320.0 * N * Math.pow (Math.cos (phi), 8.0) * l8coef * Math.pow (l, 8.0)));
234 }
235
236 /*
237 * MapXYToLatLon
238 *
239 * Converts x and y coordinates in the Transverse Mercator projection to
240 * a latitude/longitude pair. Note that Transverse Mercator is not
241 * the same as UTM; a scale factor is required to convert between them.
242 *
243 * Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J.,
244 * GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.
245 *
246 * Inputs:
247 * x - The easting of the point, in meters.
248 * y - The northing of the point, in meters.
249 * lambda0 - Longitude of the central meridian to be used, in radians.
250 *
251 * Outputs:
252 * philambda - A 2-element containing the latitude and longitude
253 * in radians.
254 *
255 * Returns:
256 * The function does not return a value.
257 *
258 * Remarks:
259 * The local variables Nf, nuf2, tf, and tf2 serve the same purpose as
260 * N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect
261 * to the footpoint latitude phif.
262 *
263 * x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and
264 * to optimize computations.
265 *
266 */
267 public LatLon MapXYToLatLon(double x, double y, double lambda0)
268 {
269 /* Get the value of phif, the footpoint latitude. */
270 double phif = FootpointLatitude (y);
271
272 /* Precalculate ep2 */
273 double ep2 = (Math.pow (Ellipsoid.GRS80.a, 2.0) - Math.pow (Ellipsoid.GRS80.b, 2.0))
274 / Math.pow (Ellipsoid.GRS80.b, 2.0);
275
276 /* Precalculate cos (phif) */
277 double cf = Math.cos (phif);
278
279 /* Precalculate nuf2 */
280 double nuf2 = ep2 * Math.pow (cf, 2.0);
281
282 /* Precalculate Nf and initialize Nfpow */
283 double Nf = Math.pow (Ellipsoid.GRS80.a, 2.0) / (Ellipsoid.GRS80.b * Math.sqrt (1 + nuf2));
284 double Nfpow = Nf;
285
286 /* Precalculate tf */
287 double tf = Math.tan (phif);
288 double tf2 = tf * tf;
289 double tf4 = tf2 * tf2;
290
291 /* Precalculate fractional coefficients for x**n in the equations
292 below to simplify the expressions for latitude and longitude. */
293 double x1frac = 1.0 / (Nfpow * cf);
294
295 Nfpow *= Nf; /* now equals Nf**2) */
296 double x2frac = tf / (2.0 * Nfpow);
297
298 Nfpow *= Nf; /* now equals Nf**3) */
299 double x3frac = 1.0 / (6.0 * Nfpow * cf);
300
301 Nfpow *= Nf; /* now equals Nf**4) */
302 double x4frac = tf / (24.0 * Nfpow);
303
304 Nfpow *= Nf; /* now equals Nf**5) */
305 double x5frac = 1.0 / (120.0 * Nfpow * cf);
306
307 Nfpow *= Nf; /* now equals Nf**6) */
308 double x6frac = tf / (720.0 * Nfpow);
309
310 Nfpow *= Nf; /* now equals Nf**7) */
311 double x7frac = 1.0 / (5040.0 * Nfpow * cf);
312
313 Nfpow *= Nf; /* now equals Nf**8) */
314 double x8frac = tf / (40320.0 * Nfpow);
315
316 /* Precalculate polynomial coefficients for x**n.
317 -- x**1 does not have a polynomial coefficient. */
318 double x2poly = -1.0 - nuf2;
319 double x3poly = -1.0 - 2 * tf2 - nuf2;
320 double x4poly = 5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2);
321 double x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2;
322 double x6poly = -61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 + 162.0 * tf2 * nuf2;
323 double x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2);
324 double x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2);
325
326 return new LatLon(
327 /* Calculate latitude */
328 Math.toDegrees(
329 phif + x2frac * x2poly * (x * x)
330 + x4frac * x4poly * Math.pow (x, 4.0)
331 + x6frac * x6poly * Math.pow (x, 6.0)
332 + x8frac * x8poly * Math.pow (x, 8.0)),
333 Math.toDegrees(
334 /* Calculate longitude */
335 lambda0 + x1frac * x
336 + x3frac * x3poly * Math.pow (x, 3.0)
337 + x5frac * x5poly * Math.pow (x, 5.0)
338 + x7frac * x7poly * Math.pow (x, 7.0)));
339 }
340
341 public EastNorth latlon2eastNorth(LatLon p) {
342 EastNorth a = MapLatLonToXY(Math.toRadians(p.lat()), Math.toRadians(p.lon()), UTMCentralMeridian(getzone()));
343 return new EastNorth(a.east() * UTMScaleFactor + 3500000.0, a.north() * UTMScaleFactor);
344 }
345
346 public LatLon eastNorth2latlon(EastNorth p) {
347 return MapXYToLatLon((p.east()-3500000.0)/UTMScaleFactor, p.north()/UTMScaleFactor, UTMCentralMeridian(getzone()));
348 }
349
350 @Override public String toString() {
351 return tr("UTM");
352 }
353
354 public int getzone()
355 {
356 return zone;
357 }
358
359 public String toCode() {
360 return "EPSG:"+ (325800 + getzone());
361 }
362
363 @Override
364 public int hashCode() {
365 return getClass().getName().hashCode()+zone; // our only real variable
366 }
367
368 public String getCacheDirectoryName() {
369 return "epsg"+ (325800 + getzone());
370 }
371
372 public Bounds getWorldBoundsLatLon()
373 {
374 return new Bounds(
375 new LatLon(-85.0, UTMCentralMeridianDeg(getzone())-5.0),
376 new LatLon(85.0, UTMCentralMeridianDeg(getzone())+5.0));
377 }
378
379 public double getDefaultZoomInPPD() {
380 // this will set the map scaler to about 1000 m
381 return 10;
382 }
383
384 public void setupPreferencePanel(JPanel p) {
385 JComboBox prefcb = new JComboBox();
386 for(int i = 1; i <= 60; i++) {
387 prefcb.addItem(i);
388 }
389
390 prefcb.setSelectedIndex(zone - 1);
391 p.setLayout(new GridBagLayout());
392 p.add(new JLabel(tr("UTM Zone")), GBC.std().insets(5,5,0,5));
393 p.add(GBC.glue(1, 0), GBC.std().fill(GBC.HORIZONTAL));
394 /* Note: we use component position 2 below to find this again */
395 p.add(prefcb, GBC.eop().fill(GBC.HORIZONTAL));
396 p.add(GBC.glue(1, 1), GBC.eol().fill(GBC.BOTH));
397 }
398
399 public Collection<String> getPreferences(JPanel p) {
400 Object prefcb = p.getComponent(2);
401 if(!(prefcb instanceof JComboBox))
402 return null;
403 int zone = ((JComboBox)prefcb).getSelectedIndex() + 1;
404 return Collections.singleton(Integer.toString(zone));
405 }
406
407 public void setPreferences(Collection<String> args)
408 {
409 zone = DEFAULT_ZONE;
410 if(args != null)
411 {
412 try {
413 for(String s : args)
414 {
415 zone = Integer.parseInt(s);
416 if(zone <= 0 || zone > 60)
417 zone = DEFAULT_ZONE;
418 break;
419 }
420 } catch(NumberFormatException e) {}
421 }
422 }
423
424 public Collection<String> getPreferencesFromCode(String code)
425 {
426 if(code.startsWith("EPSG:3258"))
427 {
428 try {
429 String zonestring = code.substring(9);
430 int zoneval = Integer.parseInt(zonestring);
431 if(zoneval > 0 && zoneval <= 60)
432 {
433 return Collections.singleton(zonestring);
434 }
435 } catch(NumberFormatException e) {}
436 }
437 return null;
438 }
439}
Note: See TracBrowser for help on using the repository browser.