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

Last change on this file since 3232 was 3083, checked in by bastiK, 14 years ago

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