001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2006, GeoTools Project Managment Committee (PMC)
005: * (C) 2001, Institut de Recherche pour le Développement
006: *
007: * This library is free software; you can redistribute it and/or
008: * modify it under the terms of the GNU Lesser General Public
009: * License as published by the Free Software Foundation;
010: * version 2.1 of the License.
011: *
012: * This library is distributed in the hope that it will be useful,
013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
015: * Lesser General Public License for more details.
016: *
017: * This class contains formulas from the public FTP area of NOAA.
018: * NOAAS's work is fully acknowledged here.
019: */
020: package org.geotools.referencing.datum;
021:
022: // J2SE dependencies and extensions
023: import java.awt.geom.Point2D;
024: import java.util.Collections;
025: import java.util.Map;
026: import javax.units.SI;
027: import javax.units.Unit;
028:
029: // OpenGIS dependencies
030: import org.opengis.referencing.datum.Ellipsoid;
031:
032: // Geotools dependencies
033: import org.geotools.geometry.GeneralDirectPosition;
034: import org.geotools.measure.CoordinateFormat;
035: import org.geotools.referencing.AbstractIdentifiedObject;
036: import org.geotools.referencing.wkt.Formatter;
037: import org.geotools.resources.Utilities;
038: import org.geotools.resources.XMath;
039: import org.geotools.resources.i18n.Errors;
040: import org.geotools.resources.i18n.ErrorKeys;
041:
042: /**
043: * Geometric figure that can be used to describe the approximate shape of the earth.
044: * In mathematical terms, it is a surface formed by the rotation of an ellipse about
045: * its minor axis. An ellipsoid requires two defining parameters:
046: * <ul>
047: * <li>{@linkplain #getSemiMajorAxis semi-major axis} and
048: * {@linkplain #getInverseFlattening inverse flattening}, or</li>
049: * <li>{@linkplain #getSemiMajorAxis semi-major axis} and
050: * {@linkplain #getSemiMinorAxis semi-minor axis}.</li>
051: * </ul>
052: *
053: * @since 2.1
054: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/datum/DefaultEllipsoid.java $
055: * @version $Id: DefaultEllipsoid.java 26412 2007-08-01 16:20:25Z aaime $
056: * @author Martin Desruisseaux
057: */
058: public class DefaultEllipsoid extends AbstractIdentifiedObject
059: implements Ellipsoid {
060: /**
061: * Serial number for interoperability with different versions.
062: */
063: private static final long serialVersionUID = -1149451543954764081L;
064:
065: /**
066: * WGS 1984 ellipsoid with axis in {@linkplain SI#METER metres}. This ellipsoid is used
067: * in GPS systems and is the default for most {@code org.geotools} packages.
068: */
069: public static final DefaultEllipsoid WGS84 = createFlattenedSphere(
070: "WGS84", 6378137.0, 298.257223563, SI.METER);
071:
072: /**
073: * GRS 80 ellipsoid with axis in {@linkplain SI#METER metres}.
074: *
075: * @since 2.2
076: */
077: public static final DefaultEllipsoid GRS80 = createFlattenedSphere(
078: "GRS80", 6378137.0, 298.257222101, SI.METER);
079:
080: /**
081: * International 1924 ellipsoid with axis in {@linkplain SI#METER metres}.
082: */
083: public static final DefaultEllipsoid INTERNATIONAL_1924 = createFlattenedSphere(
084: "International 1924", 6378388.0, 297.0, SI.METER);
085:
086: /**
087: * Clarke 1866 ellipsoid with axis in {@linkplain SI#METER metres}.
088: *
089: * @since 2.2
090: */
091: public static final DefaultEllipsoid CLARKE_1866 = createFlattenedSphere(
092: "Clarke 1866", 6378206.4, 294.9786982, SI.METER);
093:
094: /**
095: * A sphere with a radius of 6371000 {@linkplain SI#METER metres}. Spheres use a simplier
096: * algorithm for {@linkplain #orthodromicDistance orthodromic distance computation}, which
097: * may be faster and more robust.
098: */
099: public static final DefaultEllipsoid SPHERE = createEllipsoid(
100: "SPHERE", 6371000, 6371000, SI.METER);
101:
102: /**
103: * The equatorial radius.
104: * @see #getSemiMajorAxis
105: */
106: private final double semiMajorAxis;
107:
108: /**
109: * The polar radius.
110: * @see #getSemiMinorAxis
111: */
112: private final double semiMinorAxis;
113:
114: /**
115: * The inverse of the flattening value, or {@link Double#POSITIVE_INFINITY}
116: * if the ellipsoid is a sphere.
117: *
118: * @see #getInverseFlattening
119: */
120: private final double inverseFlattening;
121:
122: /**
123: * Tells if the Inverse Flattening definitive for this ellipsoid.
124: *
125: * @see #isIvfDefinitive
126: */
127: private final boolean ivfDefinitive;
128:
129: /**
130: * The units of the semi-major and semi-minor axis values.
131: */
132: private final Unit unit;
133:
134: /**
135: * Constructs a new ellipsoid with the same values than the specified one.
136: * This copy constructor provides a way to wrap an arbitrary implementation into a
137: * Geotools one or a user-defined one (as a subclass), usually in order to leverage
138: * some implementation-specific API. This constructor performs a shallow copy,
139: * i.e. the properties are not cloned.
140: *
141: * @since 2.2
142: *
143: * @see #wrap
144: */
145: protected DefaultEllipsoid(final Ellipsoid ellipsoid) {
146: super (ellipsoid);
147: semiMajorAxis = ellipsoid.getSemiMajorAxis();
148: semiMinorAxis = ellipsoid.getSemiMinorAxis();
149: inverseFlattening = ellipsoid.getInverseFlattening();
150: ivfDefinitive = ellipsoid.isIvfDefinitive();
151: unit = ellipsoid.getAxisUnit();
152: }
153:
154: /**
155: * Constructs a new ellipsoid using the specified axis length. The properties map is
156: * given unchanged to the {@linkplain AbstractIdentifiedObject#AbstractIdentifiedObject(Map)
157: * super-class constructor}.
158: *
159: * @param properties Set of properties. Should contains at least <code>"name"</code>.
160: * @param semiMajorAxis The equatorial radius.
161: * @param semiMinorAxis The polar radius.
162: * @param inverseFlattening The inverse of the flattening value.
163: * @param ivfDefinitive {@code true} if the inverse flattening is definitive.
164: * @param unit The units of the semi-major and semi-minor axis values.
165: *
166: * @see #createEllipsoid
167: * @see #createFlattenedSphere
168: */
169: protected DefaultEllipsoid(final Map properties,
170: final double semiMajorAxis, final double semiMinorAxis,
171: final double inverseFlattening,
172: final boolean ivfDefinitive, final Unit unit) {
173: super (properties);
174: this .unit = unit;
175: this .semiMajorAxis = check("semiMajorAxis", semiMajorAxis);
176: this .semiMinorAxis = check("semiMinorAxis", semiMinorAxis);
177: this .inverseFlattening = check("inverseFlattening",
178: inverseFlattening);
179: this .ivfDefinitive = ivfDefinitive;
180: ensureNonNull("unit", unit);
181: ensureLinearUnit(unit);
182: }
183:
184: /**
185: * Constructs a new ellipsoid using the specified axis length.
186: *
187: * @param name The ellipsoid name.
188: * @param semiMajorAxis The equatorial radius.
189: * @param semiMinorAxis The polar radius.
190: * @param unit The units of the semi-major and semi-minor axis values.
191: */
192: public static DefaultEllipsoid createEllipsoid(final String name,
193: final double semiMajorAxis, final double semiMinorAxis,
194: final Unit unit) {
195: return createEllipsoid(
196: Collections.singletonMap(NAME_KEY, name),
197: semiMajorAxis, semiMinorAxis, unit);
198: }
199:
200: /**
201: * Constructs a new ellipsoid using the specified axis length. The properties map is
202: * given unchanged to the {@linkplain AbstractIdentifiedObject#AbstractIdentifiedObject(Map)
203: * super-class constructor}.
204: *
205: * @param properties Set of properties. Should contains at least <code>"name"</code>.
206: * @param semiMajorAxis The equatorial radius.
207: * @param semiMinorAxis The polar radius.
208: * @param unit The units of the semi-major and semi-minor axis values.
209: */
210: public static DefaultEllipsoid createEllipsoid(
211: final Map properties, final double semiMajorAxis,
212: final double semiMinorAxis, final Unit unit) {
213: if (semiMajorAxis == semiMinorAxis) {
214: return new Spheroid(properties, semiMajorAxis, false, unit);
215: } else {
216: return new DefaultEllipsoid(properties, semiMajorAxis,
217: semiMinorAxis, semiMajorAxis
218: / (semiMajorAxis - semiMinorAxis), false,
219: unit);
220: }
221: }
222:
223: /**
224: * Constructs a new ellipsoid using the specified axis length and inverse flattening value.
225: *
226: * @param name The ellipsoid name.
227: * @param semiMajorAxis The equatorial radius.
228: * @param inverseFlattening The inverse flattening value.
229: * @param unit The units of the semi-major and semi-minor axis
230: * values.
231: */
232: public static DefaultEllipsoid createFlattenedSphere(
233: final String name, final double semiMajorAxis,
234: final double inverseFlattening, final Unit unit) {
235: return createFlattenedSphere(Collections.singletonMap(NAME_KEY,
236: name), semiMajorAxis, inverseFlattening, unit);
237: }
238:
239: /**
240: * Constructs a new ellipsoid using the specified axis length and
241: * inverse flattening value. The properties map is given unchanged to the
242: * {@linkplain AbstractIdentifiedObject#AbstractIdentifiedObject(Map) super-class constructor}.
243: *
244: * @param properties Set of properties. Should contains at least <code>"name"</code>.
245: * @param semiMajorAxis The equatorial radius.
246: * @param inverseFlattening The inverse flattening value.
247: * @param unit The units of the semi-major and semi-minor axis
248: * values.
249: */
250: public static DefaultEllipsoid createFlattenedSphere(
251: final Map properties, final double semiMajorAxis,
252: final double inverseFlattening, final Unit unit) {
253: if (Double.isInfinite(inverseFlattening)) {
254: return new Spheroid(properties, semiMajorAxis, true, unit);
255: } else {
256: return new DefaultEllipsoid(properties, semiMajorAxis,
257: semiMajorAxis * (1 - 1 / inverseFlattening),
258: inverseFlattening, true, unit);
259: }
260: }
261:
262: /**
263: * Wraps an arbitrary ellipsoid into a Geotools implementation. This method is usefull if
264: * {@link #orthodromicDistance orthodromic distance computation} (for example) are desired.
265: * If the supplied ellipsoid is already an instance of {@code DefaultEllipsoid} or is
266: * {@code null}, then it is returned unchanged.
267: */
268: public static DefaultEllipsoid wrap(final Ellipsoid ellipsoid) {
269: if (ellipsoid == null || ellipsoid instanceof DefaultEllipsoid) {
270: return (DefaultEllipsoid) ellipsoid;
271: }
272: if (ellipsoid.isIvfDefinitive()) {
273: return createFlattenedSphere(getProperties(ellipsoid),
274: ellipsoid.getSemiMajorAxis(), ellipsoid
275: .getInverseFlattening(), ellipsoid
276: .getAxisUnit());
277: } else {
278: return createEllipsoid(getProperties(ellipsoid), ellipsoid
279: .getSemiMajorAxis(), ellipsoid.getSemiMinorAxis(),
280: ellipsoid.getAxisUnit());
281: }
282: }
283:
284: /**
285: * Checks the argument validity. Argument {@code value} should be greater than zero.
286: *
287: * @param name Argument name.
288: * @param value Argument value.
289: * @return {@code value}.
290: * @throws IllegalArgumentException if {@code value} is not greater than 0.
291: */
292: static double check(final String name, final double value)
293: throws IllegalArgumentException {
294: if (value > 0) {
295: return value;
296: }
297: throw new IllegalArgumentException(Errors.format(
298: ErrorKeys.ILLEGAL_ARGUMENT_$2, name, new Double(value)));
299: }
300:
301: /**
302: * Returns the linear unit of the {@linkplain #getSemiMajorAxis semi-major}
303: * and {@linkplain #getSemiMinorAxis semi-minor} axis values.
304: *
305: * @return The axis linear unit.
306: */
307: public Unit getAxisUnit() {
308: return unit;
309: }
310:
311: /**
312: * Length of the semi-major axis of the ellipsoid. This is the
313: * equatorial radius in {@linkplain #getAxisUnit axis linear unit}.
314: *
315: * @return Length of semi-major axis.
316: */
317: public double getSemiMajorAxis() {
318: return semiMajorAxis;
319: }
320:
321: /**
322: * Length of the semi-minor axis of the ellipsoid. This is the
323: * polar radius in {@linkplain #getAxisUnit axis linear unit}.
324: *
325: * @return Length of semi-minor axis.
326: */
327: public double getSemiMinorAxis() {
328: return semiMinorAxis;
329: }
330:
331: /**
332: * The ratio of the distance between the center and a focus of the ellipse
333: * to the length of its semimajor axis. The eccentricity can alternately be
334: * computed from the equation: <code>e=sqrt(2f-f˛)</code>.
335: */
336: public double getEccentricity() {
337: final double f = 1 - getSemiMinorAxis() / getSemiMajorAxis();
338: return Math.sqrt(2 * f - f * f);
339: }
340:
341: /**
342: * Returns the value of the inverse of the flattening constant. Flattening is a value
343: * used to indicate how closely an ellipsoid approaches a spherical shape. The inverse
344: * flattening is related to the equatorial/polar radius by the formula
345: *
346: * <var>ivf</var> = <var>r</var><sub>e</sub>/(<var>r</var><sub>e</sub>-<var>r</var><sub>p</sub>).
347: *
348: * For perfect spheres (i.e. if {@link #isSphere} returns {@code true}),
349: * the {@link Double#POSITIVE_INFINITY} value is used.
350: *
351: * @return The inverse flattening value.
352: */
353: public double getInverseFlattening() {
354: return inverseFlattening;
355: }
356:
357: /**
358: * Indicates if the {@linkplain #getInverseFlattening inverse flattening} is definitive for
359: * this ellipsoid. Some ellipsoids use the IVF as the defining value, and calculate the polar
360: * radius whenever asked. Other ellipsoids use the polar radius to calculate the IVF whenever
361: * asked. This distinction can be important to avoid floating-point rounding errors.
362: *
363: * @return {@code true} if the {@linkplain #getInverseFlattening inverse flattening} is
364: * definitive, or {@code false} if the {@linkplain #getSemiMinorAxis polar radius}
365: * is definitive.
366: */
367: public boolean isIvfDefinitive() {
368: return ivfDefinitive;
369: }
370:
371: /**
372: * {@code true} if the ellipsoid is degenerate and is actually a sphere. The sphere is
373: * completely defined by the {@linkplain #getSemiMajorAxis semi-major axis}, which is the
374: * radius of the sphere.
375: *
376: * @return {@code true} if the ellipsoid is degenerate and is actually a sphere.
377: */
378: public boolean isSphere() {
379: return semiMajorAxis == semiMinorAxis;
380: }
381:
382: /**
383: * Returns the orthodromic distance between two geographic coordinates.
384: * The orthodromic distance is the shortest distance between two points
385: * on a sphere's surface. The default implementation delegates the work
386: * to {@link #orthodromicDistance(double,double,double,double)}.
387: *
388: * @param P1 Longitude and latitude of first point (in decimal degrees).
389: * @param P2 Longitude and latitude of second point (in decimal degrees).
390: * @return The orthodromic distance (in the units of this ellipsoid).
391: */
392: public double orthodromicDistance(final Point2D P1, final Point2D P2) {
393: return orthodromicDistance(P1.getX(), P1.getY(), P2.getX(), P2
394: .getY());
395: }
396:
397: /**
398: * Returns the orthodromic distance between two geographic coordinates.
399: * The orthodromic distance is the shortest distance between two points
400: * on a sphere's surface. The orthodromic path is always on a great circle.
401: * This is different from the <cite>loxodromic distance</cite>, which is a
402: * longer distance on a path with a constant direction on the compass.
403: *
404: * @param x1 Longitude of first point (in decimal degrees).
405: * @param y1 Latitude of first point (in decimal degrees).
406: * @param x2 Longitude of second point (in decimal degrees).
407: * @param y2 Latitude of second point (in decimal degrees).
408: * @return The orthodromic distance (in the units of this ellipsoid's axis).
409: */
410: public double orthodromicDistance(double x1, double y1, double x2,
411: double y2) {
412: x1 = Math.toRadians(x1);
413: y1 = Math.toRadians(y1);
414: x2 = Math.toRadians(x2);
415: y2 = Math.toRadians(y2);
416: /*
417: * Solution of the geodetic inverse problem after T.Vincenty.
418: * Modified Rainsford's method with Helmert's elliptical terms.
419: * Effective in any azimuth and at any distance short of antipodal.
420: *
421: * Latitudes and longitudes in radians positive North and East.
422: * Forward azimuths at both points returned in radians from North.
423: *
424: * Programmed for CDC-6600 by LCDR L.Pfeifer NGS ROCKVILLE MD 18FEB75
425: * Modified for IBM SYSTEM 360 by John G.Gergen NGS ROCKVILLE MD 7507
426: * Ported from Fortran to Java by Martin Desruisseaux.
427: *
428: * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/inverse.for
429: * subroutine INVER1
430: */
431: final int MAX_ITERATIONS = 100;
432: final double EPS = 0.5E-13;
433: final double F = 1 / getInverseFlattening();
434: final double R = 1 - F;
435:
436: double tu1 = R * Math.sin(y1) / Math.cos(y1);
437: double tu2 = R * Math.sin(y2) / Math.cos(y2);
438: double cu1 = 1 / Math.sqrt(tu1 * tu1 + 1);
439: double cu2 = 1 / Math.sqrt(tu2 * tu2 + 1);
440: double su1 = cu1 * tu1;
441: double s = cu1 * cu2;
442: double baz = s * tu2;
443: double faz = baz * tu1;
444: double x = x2 - x1;
445: for (int i = 0; i < MAX_ITERATIONS; i++) {
446: final double sx = Math.sin(x);
447: final double cx = Math.cos(x);
448: tu1 = cu2 * sx;
449: tu2 = baz - su1 * cu2 * cx;
450: final double sy = XMath.hypot(tu1, tu2);
451: final double cy = s * cx + faz;
452: final double y = Math.atan2(sy, cy);
453: final double SA = s * sx / sy;
454: final double c2a = 1 - SA * SA;
455: double cz = faz + faz;
456: if (c2a > 0) {
457: cz = -cz / c2a + cy;
458: }
459: double e = cz * cz * 2 - 1;
460: double c = ((-3 * c2a + 4) * F + 4) * c2a * F / 16;
461: double d = x;
462: x = ((e * cy * c + cz) * sy * c + y) * SA;
463: x = (1 - c) * x * F + x2 - x1;
464:
465: if (Math.abs(d - x) <= EPS) {
466: if (false) {
467: // 'faz' and 'baz' are forward azimuths at both points.
468: // Since the current API can't returns this result, it
469: // doesn't worth to compute it at this time.
470: faz = Math.atan2(tu1, tu2);
471: baz = Math.atan2(cu1 * sx, baz * cx - su1 * cu2)
472: + Math.PI;
473: }
474: x = Math.sqrt((1 / (R * R) - 1) * c2a + 1) + 1;
475: x = (x - 2) / x;
476: c = 1 - x;
477: c = (x * x / 4 + 1) / c;
478: d = (0.375 * x * x - 1) * x;
479: x = e * cy;
480: s = 1 - 2 * e;
481: s = ((((sy * sy * 4 - 3) * s * cz * d / 6 - x) * d / 4 + cz)
482: * sy * d + y)
483: * c * R * getSemiMajorAxis();
484: return s;
485: }
486: }
487: // No convergence. It may be because coordinate points
488: // are equals or because they are at antipodes.
489: final double LEPS = 1E-10;
490: if (Math.abs(x1 - x2) <= LEPS && Math.abs(y1 - y2) <= LEPS) {
491: return 0; // Coordinate points are equals
492: }
493: if (Math.abs(y1) <= LEPS && Math.abs(y2) <= LEPS) {
494: return Math.abs(x1 - x2) * getSemiMajorAxis(); // Points are on the equator.
495: }
496: // Other cases: no solution for this algorithm.
497: final CoordinateFormat format = new CoordinateFormat();
498: throw new ArithmeticException(Errors.format(
499: ErrorKeys.NO_CONVERGENCE_$2, format
500: .format(new GeneralDirectPosition(Math
501: .toDegrees(x1), Math.toDegrees(y1))),
502: format.format(new GeneralDirectPosition(Math
503: .toDegrees(x2), Math.toDegrees(y2)))));
504: }
505:
506: /**
507: * Compare this ellipsoid with the specified object for equality.
508: *
509: * @param object The object to compare to {@code this}.
510: * @param compareMetadata {@code true} for performing a strict comparaison, or
511: * {@code false} for comparing only properties relevant to transformations.
512: * @return {@code true} if both objects are equal.
513: */
514: public boolean equals(final AbstractIdentifiedObject object,
515: final boolean compareMetadata) {
516: if (object == this ) {
517: return true; // Slight optimization.
518: }
519: if (super .equals(object, compareMetadata)) {
520: final DefaultEllipsoid that = (DefaultEllipsoid) object;
521: return (!compareMetadata || this .ivfDefinitive == that.ivfDefinitive)
522: && Double.doubleToLongBits(this .semiMajorAxis) == Double
523: .doubleToLongBits(that.semiMajorAxis)
524: && Double.doubleToLongBits(this .semiMinorAxis) == Double
525: .doubleToLongBits(that.semiMinorAxis)
526: && Double.doubleToLongBits(this .inverseFlattening) == Double
527: .doubleToLongBits(that.inverseFlattening)
528: && Utilities.equals(this .unit, that.unit);
529: }
530: return false;
531: }
532:
533: /**
534: * Returns a hash value for this ellipsoid. {@linkplain #getName Name},
535: * {@linkplain #getRemarks remarks} and the like are not taken in account.
536: * In other words, two ellipsoids will return the same hash value if they
537: * are equal in the sense of
538: * <code>{@link #equals equals}(AbstractIdentifiedObject, <strong>false</strong>)</code>.
539: *
540: * @return The hash code value. This value doesn't need to be the same
541: * in past or future versions of this class.
542: */
543: public int hashCode() {
544: long longCode = 37 * Double.doubleToLongBits(semiMajorAxis);
545: if (ivfDefinitive) {
546: longCode += inverseFlattening;
547: } else {
548: longCode += semiMinorAxis;
549: }
550: return (((int) (longCode >>> 32)) ^ (int) longCode);
551: }
552:
553: /**
554: * Format the inner part of a
555: * <A HREF="http://geoapi.sourceforge.net/snapshot/javadoc/org/opengis/referencing/doc-files/WKT.html"><cite>Well
556: * Known Text</cite> (WKT)</A> element.
557: *
558: * @param formatter The formatter to use.
559: * @return The WKT element name, which is "SPHEROID"
560: */
561: protected String formatWKT(final Formatter formatter) {
562: final double ivf = getInverseFlattening();
563: formatter.append(getAxisUnit().getConverterTo(SI.METER)
564: .convert(getSemiMajorAxis()));
565: formatter.append(Double.isInfinite(ivf) ? 0 : ivf);
566: return "SPHEROID";
567: }
568: }
|