1 | /*
|
---|
2 | * Import from fr.geo.convert package, a geographic coordinates converter.
|
---|
3 | * (https://www.i3s.unice.fr/~johan/gps/)
|
---|
4 | * License: GPL. For details, see LICENSE file.
|
---|
5 | * Copyright (C) 2002 Johan Montagnat (johan@creatis.insa-lyon.fr)
|
---|
6 | */
|
---|
7 | package org.openstreetmap.josm.data.projection;
|
---|
8 |
|
---|
9 | import org.openstreetmap.josm.data.coor.LatLon;
|
---|
10 | import org.openstreetmap.josm.tools.Utils;
|
---|
11 |
|
---|
12 | /**
|
---|
13 | * Reference ellipsoids.
|
---|
14 | */
|
---|
15 | public final class Ellipsoid {
|
---|
16 |
|
---|
17 | /**
|
---|
18 | * Airy 1830
|
---|
19 | */
|
---|
20 | public static final Ellipsoid Airy = Ellipsoid.createAb(6377563.396, 6356256.910);
|
---|
21 |
|
---|
22 | /**
|
---|
23 | * Modified Airy 1849
|
---|
24 | */
|
---|
25 | public static final Ellipsoid AiryMod = Ellipsoid.createAb(6377340.189, 6356034.446);
|
---|
26 |
|
---|
27 | /**
|
---|
28 | * Australian National Spheroid (Australian Natl & S. Amer. 1969)
|
---|
29 | * same as GRS67 Modified
|
---|
30 | */
|
---|
31 | public static final Ellipsoid AustSA = Ellipsoid.createArf(6378160.0, 298.25);
|
---|
32 |
|
---|
33 | /**
|
---|
34 | * Bessel 1841 ellipsoid
|
---|
35 | */
|
---|
36 | public static final Ellipsoid Bessel1841 = Ellipsoid.createArf(6377397.155, 299.1528128);
|
---|
37 |
|
---|
38 | /**
|
---|
39 | * Bessel 1841 (Namibia)
|
---|
40 | */
|
---|
41 | public static final Ellipsoid BesselNamibia = Ellipsoid.createArf(6377483.865, 299.1528128);
|
---|
42 |
|
---|
43 | /**
|
---|
44 | * Clarke 1866 ellipsoid
|
---|
45 | */
|
---|
46 | public static final Ellipsoid Clarke1866 = Ellipsoid.createAb(6378206.4, 6356583.8);
|
---|
47 |
|
---|
48 | /**
|
---|
49 | * Clarke 1880 (modified)
|
---|
50 | */
|
---|
51 | public static final Ellipsoid Clarke1880 = Ellipsoid.createArf(6378249.145, 293.4663);
|
---|
52 |
|
---|
53 | /**
|
---|
54 | * Clarke 1880 IGN (French national geographic institute)
|
---|
55 | */
|
---|
56 | public static final Ellipsoid ClarkeIGN = Ellipsoid.createAb(6378249.2, 6356515.0);
|
---|
57 |
|
---|
58 | /**
|
---|
59 | * Everest 1830
|
---|
60 | */
|
---|
61 | public static final Ellipsoid Everest = Ellipsoid.createArf(6377276.345, 300.8017);
|
---|
62 |
|
---|
63 | /**
|
---|
64 | * Everest 1948
|
---|
65 | */
|
---|
66 | public static final Ellipsoid Everest1948 = Ellipsoid.createArf(6377304.063, 300.8017);
|
---|
67 |
|
---|
68 | /**
|
---|
69 | * Everest 1956
|
---|
70 | */
|
---|
71 | public static final Ellipsoid Everest1956 = Ellipsoid.createArf(6377301.243, 300.8017);
|
---|
72 |
|
---|
73 | /**
|
---|
74 | * Everest 1969
|
---|
75 | */
|
---|
76 | public static final Ellipsoid Everest1969 = Ellipsoid.createArf(6377295.664, 300.8017);
|
---|
77 |
|
---|
78 | /**
|
---|
79 | * Everest (Sabah & Sarawak)
|
---|
80 | */
|
---|
81 | public static final Ellipsoid EverestSabahSarawak = Ellipsoid.createArf(6377298.556, 300.8017);
|
---|
82 |
|
---|
83 | /**
|
---|
84 | * Fischer (Mercury Datum) 1960
|
---|
85 | */
|
---|
86 | public static final Ellipsoid Fischer = Ellipsoid.createArf(6378166., 298.3);
|
---|
87 |
|
---|
88 | /**
|
---|
89 | * Modified Fischer 1960
|
---|
90 | */
|
---|
91 | public static final Ellipsoid FischerMod = Ellipsoid.createArf(6378155., 298.3);
|
---|
92 |
|
---|
93 | /**
|
---|
94 | * Fischer 1968
|
---|
95 | */
|
---|
96 | public static final Ellipsoid Fischer1968 = Ellipsoid.createArf(6378150., 298.3);
|
---|
97 |
|
---|
98 | /**
|
---|
99 | * GRS67 ellipsoid
|
---|
100 | */
|
---|
101 | public static final Ellipsoid GRS67 = Ellipsoid.createArf(6378160.0, 298.247167427);
|
---|
102 |
|
---|
103 | /**
|
---|
104 | * GRS80 ellipsoid
|
---|
105 | */
|
---|
106 | public static final Ellipsoid GRS80 = Ellipsoid.createArf(6378137.0, 298.257222101);
|
---|
107 |
|
---|
108 | /**
|
---|
109 | * Hayford's ellipsoid 1909 (ED50 system)
|
---|
110 | * Also known as International 1924
|
---|
111 | * Proj.4 code: intl
|
---|
112 | */
|
---|
113 | public static final Ellipsoid Hayford = Ellipsoid.createArf(6378388.0, 297.0);
|
---|
114 |
|
---|
115 | /**
|
---|
116 | * Helmert 1906
|
---|
117 | */
|
---|
118 | public static final Ellipsoid Helmert = Ellipsoid.createArf(6378200.0, 298.3);
|
---|
119 |
|
---|
120 | /**
|
---|
121 | * Hough
|
---|
122 | */
|
---|
123 | public static final Ellipsoid Hough = Ellipsoid.createArf(6378270.0, 297.0);
|
---|
124 |
|
---|
125 | /**
|
---|
126 | * Krassowsky 1940 ellipsoid
|
---|
127 | */
|
---|
128 | public static final Ellipsoid Krassowsky = Ellipsoid.createArf(6378245.0, 298.3);
|
---|
129 |
|
---|
130 | /**
|
---|
131 | * Sphere
|
---|
132 | */
|
---|
133 | public static final Ellipsoid Sphere = Ellipsoid.createAb(6370997.0, 6370997.0);
|
---|
134 |
|
---|
135 | /**
|
---|
136 | * Walbeck
|
---|
137 | */
|
---|
138 | public static final Ellipsoid Walbeck = Ellipsoid.createAb(6376896.0, 6355834.8467);
|
---|
139 |
|
---|
140 | /**
|
---|
141 | * WGS66 ellipsoid
|
---|
142 | */
|
---|
143 | public static final Ellipsoid WGS66 = Ellipsoid.createArf(6378145.0, 298.25);
|
---|
144 |
|
---|
145 | /**
|
---|
146 | * WGS72 ellipsoid
|
---|
147 | */
|
---|
148 | public static final Ellipsoid WGS72 = Ellipsoid.createArf(6378135.0, 298.26);
|
---|
149 |
|
---|
150 | /**
|
---|
151 | * WGS84 ellipsoid
|
---|
152 | */
|
---|
153 | public static final Ellipsoid WGS84 = Ellipsoid.createArf(6378137.0, 298.257223563);
|
---|
154 |
|
---|
155 | /**
|
---|
156 | * half long axis
|
---|
157 | */
|
---|
158 | public final double a;
|
---|
159 |
|
---|
160 | /**
|
---|
161 | * half short axis
|
---|
162 | */
|
---|
163 | public final double b;
|
---|
164 |
|
---|
165 | /**
|
---|
166 | * first eccentricity:
|
---|
167 | * sqrt(a*a - b*b) / a
|
---|
168 | */
|
---|
169 | public final double e;
|
---|
170 |
|
---|
171 | /**
|
---|
172 | * first eccentricity squared:
|
---|
173 | * (a*a - b*b) / (a*a)
|
---|
174 | */
|
---|
175 | public final double e2;
|
---|
176 |
|
---|
177 | /**
|
---|
178 | * square of the second eccentricity:
|
---|
179 | * (a*a - b*b) / (b*b)
|
---|
180 | */
|
---|
181 | public final double eb2;
|
---|
182 |
|
---|
183 | /**
|
---|
184 | * if ellipsoid is spherical, i.e. the major and minor semiaxis are
|
---|
185 | * the same
|
---|
186 | */
|
---|
187 | public final boolean spherical;
|
---|
188 |
|
---|
189 | /**
|
---|
190 | * private constructur - use one of the create_* methods
|
---|
191 | *
|
---|
192 | * @param a semimajor radius of the ellipsoid axis
|
---|
193 | * @param b semiminor radius of the ellipsoid axis
|
---|
194 | * @param e first eccentricity of the ellipsoid ( = sqrt((a*a - b*b)/(a*a)))
|
---|
195 | * @param e2 first eccentricity squared
|
---|
196 | * @param eb2 square of the second eccentricity
|
---|
197 | * @param sperical if the ellipsoid is sphere
|
---|
198 | */
|
---|
199 | private Ellipsoid(double a, double b, double e, double e2, double eb2, boolean sperical) {
|
---|
200 | this.a = a;
|
---|
201 | this.b = b;
|
---|
202 | this.e = e;
|
---|
203 | this.e2 = e2;
|
---|
204 | this.eb2 = eb2;
|
---|
205 | this.spherical = sperical;
|
---|
206 | }
|
---|
207 |
|
---|
208 | /**
|
---|
209 | * create a new ellipsoid
|
---|
210 | *
|
---|
211 | * @param a semimajor radius of the ellipsoid axis (in meters)
|
---|
212 | * @param b semiminor radius of the ellipsoid axis (in meters)
|
---|
213 | * @return the new ellipsoid
|
---|
214 | */
|
---|
215 | public static Ellipsoid createAb(double a, double b) {
|
---|
216 | double e2 = (a*a - b*b) / (a*a);
|
---|
217 | double e = Math.sqrt(e2);
|
---|
218 | double eb2 = e2 / (1.0 - e2);
|
---|
219 | return new Ellipsoid(a, b, e, e2, eb2, a == b);
|
---|
220 | }
|
---|
221 |
|
---|
222 | /**
|
---|
223 | * create a new ellipsoid
|
---|
224 | *
|
---|
225 | * @param a semimajor radius of the ellipsoid axis (in meters)
|
---|
226 | * @param es first eccentricity squared
|
---|
227 | * @return the new ellipsoid
|
---|
228 | */
|
---|
229 | public static Ellipsoid createAes(double a, double es) {
|
---|
230 | double b = a * Math.sqrt(1.0 - es);
|
---|
231 | double e = Math.sqrt(es);
|
---|
232 | double eb2 = es / (1.0 - es);
|
---|
233 | return new Ellipsoid(a, b, e, es, eb2, es == 0);
|
---|
234 | }
|
---|
235 |
|
---|
236 | /**
|
---|
237 | * create a new ellipsoid
|
---|
238 | *
|
---|
239 | * @param a semimajor radius of the ellipsoid axis (in meters)
|
---|
240 | * @param f flattening ( = (a - b) / a)
|
---|
241 | * @return the new ellipsoid
|
---|
242 | */
|
---|
243 | public static Ellipsoid createAf(double a, double f) {
|
---|
244 | double b = a * (1.0 - f);
|
---|
245 | double e2 = f * (2 - f);
|
---|
246 | double e = Math.sqrt(e2);
|
---|
247 | double eb2 = e2 / (1.0 - e2);
|
---|
248 | return new Ellipsoid(a, b, e, e2, eb2, f == 0);
|
---|
249 | }
|
---|
250 |
|
---|
251 | /**
|
---|
252 | * create a new ellipsoid
|
---|
253 | *
|
---|
254 | * @param a semimajor radius of the ellipsoid axis (in meters)
|
---|
255 | * @param rf inverse flattening
|
---|
256 | * @return the new ellipsoid
|
---|
257 | */
|
---|
258 | public static Ellipsoid createArf(double a, double rf) {
|
---|
259 | return createAf(a, 1.0 / rf);
|
---|
260 | }
|
---|
261 |
|
---|
262 | @Override
|
---|
263 | public String toString() {
|
---|
264 | return "Ellipsoid{a="+a+", b="+b+'}';
|
---|
265 | }
|
---|
266 |
|
---|
267 | /**
|
---|
268 | * Returns the <i>radius of curvature in the prime vertical</i>
|
---|
269 | * for this reference ellipsoid at the specified latitude.
|
---|
270 | *
|
---|
271 | * @param phi The local latitude (radians).
|
---|
272 | * @return The radius of curvature in the prime vertical (meters).
|
---|
273 | */
|
---|
274 | public double verticalRadiusOfCurvature(final double phi) {
|
---|
275 | return a / Math.sqrt(1.0 - (e2 * sqr(Math.sin(phi))));
|
---|
276 | }
|
---|
277 |
|
---|
278 | private static double sqr(final double x) {
|
---|
279 | return x * x;
|
---|
280 | }
|
---|
281 |
|
---|
282 | /**
|
---|
283 | * Returns the meridional arc, the true meridional distance on the
|
---|
284 | * ellipsoid from the equator to the specified latitude, in meters.
|
---|
285 | *
|
---|
286 | * @param phi The local latitude (in radians).
|
---|
287 | * @return The meridional arc (in meters).
|
---|
288 | */
|
---|
289 | public double meridionalArc(final double phi) {
|
---|
290 | final double sin2Phi = Math.sin(2.0 * phi);
|
---|
291 | final double sin4Phi = Math.sin(4.0 * phi);
|
---|
292 | final double sin6Phi = Math.sin(6.0 * phi);
|
---|
293 | final double sin8Phi = Math.sin(8.0 * phi);
|
---|
294 | // TODO . calculate 'f'
|
---|
295 | //double f = 1.0 / 298.257222101; // GRS80
|
---|
296 | double f = 1.0 / 298.257223563; // WGS84
|
---|
297 | final double n = f / (2.0 - f);
|
---|
298 | final double n2 = n * n;
|
---|
299 | final double n3 = n2 * n;
|
---|
300 | final double n4 = n3 * n;
|
---|
301 | final double n5 = n4 * n;
|
---|
302 | final double n1n2 = n - n2;
|
---|
303 | final double n2n3 = n2 - n3;
|
---|
304 | final double n3n4 = n3 - n4;
|
---|
305 | final double n4n5 = n4 - n5;
|
---|
306 | final double ap = a * (1.0 - n + (5.0 / 4.0) * n2n3 + (81.0 / 64.0) * n4n5);
|
---|
307 | final double bp = (3.0 / 2.0) * a * (n1n2 + (7.0 / 8.0) * n3n4 + (55.0 / 64.0) * n5);
|
---|
308 | final double cp = (15.0 / 16.0) * a * (n2n3 + (3.0 / 4.0) * n4n5);
|
---|
309 | final double dp = (35.0 / 48.0) * a * (n3n4 + (11.0 / 16.0) * n5);
|
---|
310 | final double ep = (315.0 / 512.0) * a * n4n5;
|
---|
311 | return ap * phi - bp * sin2Phi + cp * sin4Phi - dp * sin6Phi + ep * sin8Phi;
|
---|
312 | }
|
---|
313 |
|
---|
314 | /**
|
---|
315 | * Returns the <i>radius of curvature in the meridian</i>
|
---|
316 | * for this reference ellipsoid at the specified latitude.
|
---|
317 | *
|
---|
318 | * @param phi The local latitude (in radians).
|
---|
319 | * @return The radius of curvature in the meridian (in meters).
|
---|
320 | */
|
---|
321 | public double meridionalRadiusOfCurvature(final double phi) {
|
---|
322 | return verticalRadiusOfCurvature(phi)
|
---|
323 | / (1.0 + eb2 * sqr(Math.cos(phi)));
|
---|
324 | }
|
---|
325 |
|
---|
326 | /**
|
---|
327 | * Returns isometric latitude of phi on given first eccentricity (e)
|
---|
328 | * @param phi The local latitude (radians).
|
---|
329 | * @param e first eccentricity
|
---|
330 | * @return isometric latitude of phi on first eccentricity (e)
|
---|
331 | */
|
---|
332 | public double latitudeIsometric(double phi, double e) {
|
---|
333 | double v1 = 1-e*Math.sin(phi);
|
---|
334 | double v2 = 1+e*Math.sin(phi);
|
---|
335 | return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2));
|
---|
336 | }
|
---|
337 |
|
---|
338 | /**
|
---|
339 | * Returns isometric latitude of phi on first eccentricity (e)
|
---|
340 | * @param phi The local latitude (radians).
|
---|
341 | * @return isometric latitude of phi on first eccentricity (e)
|
---|
342 | */
|
---|
343 | public double latitudeIsometric(double phi) {
|
---|
344 | double v1 = 1-e*Math.sin(phi);
|
---|
345 | double v2 = 1+e*Math.sin(phi);
|
---|
346 | return Math.log(Math.tan(Math.PI/4+phi/2)*Math.pow(v1/v2, e/2));
|
---|
347 | }
|
---|
348 |
|
---|
349 | /**
|
---|
350 | * Returns geographic latitude of isometric latitude of first eccentricity (e) and epsilon precision
|
---|
351 | * @param latIso isometric latitude
|
---|
352 | * @param e first eccentricity
|
---|
353 | * @param epsilon epsilon precision
|
---|
354 | * @return geographic latitude of isometric latitude of first eccentricity (e) and epsilon precision
|
---|
355 | */
|
---|
356 | public double latitude(double latIso, double e, double epsilon) {
|
---|
357 | double lat0 = 2*Math.atan(Math.exp(latIso))-Math.PI/2;
|
---|
358 | double lati = lat0;
|
---|
359 | double lati1 = 1.0; // random value to start the iterative processes
|
---|
360 | while (Math.abs(lati1-lati) >= epsilon) {
|
---|
361 | lati = lati1;
|
---|
362 | double v1 = 1+e*Math.sin(lati);
|
---|
363 | double v2 = 1-e*Math.sin(lati);
|
---|
364 | lati1 = 2*Math.atan(Math.pow(v1/v2, e/2)*Math.exp(latIso))-Math.PI/2;
|
---|
365 | }
|
---|
366 | return lati1;
|
---|
367 | }
|
---|
368 |
|
---|
369 | /**
|
---|
370 | * convert cartesian coordinates to ellipsoidal coordinates
|
---|
371 | *
|
---|
372 | * @param xyz the coordinates in meters (X, Y, Z)
|
---|
373 | * @return The corresponding latitude and longitude in degrees
|
---|
374 | */
|
---|
375 | public LatLon cart2LatLon(double... xyz) {
|
---|
376 | return cart2LatLon(xyz, 1e-11);
|
---|
377 | }
|
---|
378 |
|
---|
379 | public LatLon cart2LatLon(double[] xyz, double epsilon) {
|
---|
380 | double norm = Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
|
---|
381 | double lg = 2.0 * Math.atan(xyz[1] / (xyz[0] + norm));
|
---|
382 | double lt = Math.atan(xyz[2] / (norm * (1.0 - (a * e2 / Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2])))));
|
---|
383 | double delta = 1.0;
|
---|
384 | while (delta > epsilon) {
|
---|
385 | double s2 = Math.sin(lt);
|
---|
386 | s2 *= s2;
|
---|
387 | double l = Math.atan((xyz[2] / norm)
|
---|
388 | / (1.0 - (a * e2 * Math.cos(lt) / (norm * Math.sqrt(1.0 - e2 * s2)))));
|
---|
389 | delta = Math.abs(l - lt);
|
---|
390 | lt = l;
|
---|
391 | }
|
---|
392 | return new LatLon(Utils.toDegrees(lt), Utils.toDegrees(lg));
|
---|
393 | }
|
---|
394 |
|
---|
395 | /**
|
---|
396 | * convert ellipsoidal coordinates to cartesian coordinates
|
---|
397 | *
|
---|
398 | * @param coord The Latitude and longitude in degrees
|
---|
399 | * @return the corresponding (X, Y Z) cartesian coordinates in meters.
|
---|
400 | */
|
---|
401 | public double[] latLon2Cart(LatLon coord) {
|
---|
402 | double phi = Utils.toRadians(coord.lat());
|
---|
403 | double lambda = Utils.toRadians(coord.lon());
|
---|
404 |
|
---|
405 | double rn = a / Math.sqrt(1 - e2 * Math.pow(Math.sin(phi), 2));
|
---|
406 | double[] xyz = new double[3];
|
---|
407 | xyz[0] = rn * Math.cos(phi) * Math.cos(lambda);
|
---|
408 | xyz[1] = rn * Math.cos(phi) * Math.sin(lambda);
|
---|
409 | xyz[2] = rn * (1 - e2) * Math.sin(phi);
|
---|
410 |
|
---|
411 | return xyz;
|
---|
412 | }
|
---|
413 | }
|
---|