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

Last change on this file since 12163 was 12161, checked in by michael2402, 7 years ago

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