source: josm/trunk/src/org/openstreetmap/josm/data/projection/TransverseMercator.java@ 3779

Last change on this file since 3779 was 3635, checked in by bastiK, 13 years ago

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