source: josm/trunk/src/org/openstreetmap/josm/data/projection/Ellipsoid.java@ 14159

Last change on this file since 14159 was 13626, checked in by Don-vip, 6 years ago

see #16129 - fix Sphere ellipsoid radius (values of proj.4 and Proj4J are different!)
  • Property svn:eol-style set to native
File size: 12.6 KB
Line 
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 */
7package org.openstreetmap.josm.data.projection;
8
9import org.openstreetmap.josm.data.coor.LatLon;
10import org.openstreetmap.josm.tools.Utils;
11
12/**
13 * Reference ellipsoids.
14 */
15public 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 processus
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}
Note: See TracBrowser for help on using the repository browser.