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

Last change on this file since 4382 was 4285, checked in by bastiK, 13 years ago

major projection rework

More modular structure, inspired by Proj.4.

There are almost no semantic changes to the projection algorithms. Mostly factors of 'a' and 180/PI have been moved from one place to the other. In UTM_France_DOM, the ellipsoid conversion for the first 3 projections has been changed from hayford <-> GRS80 to hayford <-> WGS84.

Some redundant algorithms have been removed. In particular:

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