001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2006, Geotools Project Managment Committee (PMC)
005: * (C) 2002, 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; either
010: * version 2.1 of the License, or (at your option) any later version.
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: package org.geotools.referencing.operation.transform;
018:
019: // J2SE dependencies and extensions
020: import javax.vecmath.Point3d;
021:
022: // JUnit dependencies
023: import junit.framework.Test;
024: import junit.framework.TestSuite;
025:
026: // OpenGIS dependencies
027: import org.opengis.referencing.FactoryException;
028: import org.opengis.referencing.crs.CoordinateReferenceSystem;
029: import org.opengis.referencing.operation.CoordinateOperation;
030: import org.opengis.referencing.operation.MathTransform;
031: import org.opengis.referencing.operation.TransformException;
032:
033: // Geotools dependencies
034: import org.geotools.referencing.crs.DefaultGeocentricCRS;
035: import org.geotools.referencing.crs.DefaultGeographicCRS;
036: import org.geotools.referencing.datum.DefaultEllipsoid;
037: import org.geotools.referencing.operation.TestTransform;
038: import org.geotools.resources.XMath;
039:
040: /**
041: * Test the following transformation classes with the geocentric transform:
042: *
043: * <ul>
044: * <li>{@link CoordinateOperation}</li>
045: * <li>{@link GeocentricTransform}</li>
046: * <li>{@link DefaultEllipsoid}</li>
047: * </ul>
048: *
049: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/test/java/org/geotools/referencing/operation/transform/GeocentricTransformTest.java $
050: * @version $Id: GeocentricTransformTest.java 25262 2007-04-23 21:11:16Z desruisseaux $
051: * @author Martin Desruisseaux
052: */
053: public final class GeocentricTransformTest extends TestTransform {
054: /**
055: * Runs the tests with the textual test runner.
056: */
057: public static void main(String args[]) {
058: junit.textui.TestRunner.run(suite());
059: }
060:
061: /**
062: * Returns the test suite.
063: */
064: public static Test suite() {
065: return new TestSuite(GeocentricTransformTest.class);
066: }
067:
068: /**
069: * Constructs a test case with the given name.
070: */
071: public GeocentricTransformTest(final String name) {
072: super (name);
073: }
074:
075: /**
076: * Tests the orthodromic distance computed by {@link DefaultEllipsoid}. There is actually two
077: * algorithms used: one for the ellipsoidal model, and a simpler one for spherical model.
078: * We test the ellipsoidal model using know values of nautical mile at different latitude.
079: * Then, we test the spherical model with random values. If JDK 1.4 assertion is enabled,
080: * the spherical model will compare its result with the ellipsoidal one.
081: *
082: * Note about nautical mile:
083: *
084: * "Le mille marin était, en principe, la longueur de la minute sexagésimale du méridien
085: * à la latitude de 45°. Cette longueur dépendait donc des valeurs adoptées pour le rayon
086: * équatorial de la terre et son aplatissement. En France, le décret du 3 mai 1961 sur les
087: * unités de mesure, fixe à 1852 mètres la longueur du mille marin qui est également la
088: * valeur adoptée pour le mille marin international."
089: *
090: * Source: Office de la langue française, 1996
091: * http://www.granddictionnaire.com
092: */
093: public void testEllipsoid() throws FactoryException {
094: final DefaultEllipsoid e = DefaultEllipsoid.WGS84;
095: final double hm = 0.5 / 60; // Half of a minute of angle, in degrees.
096: /*
097: * Test the ellipsoidal model.
098: */
099: assertEquals("Nautical mile at equator", 1842.78, e
100: .orthodromicDistance(0, 00 - hm, 0, 00 + hm), 0.2);
101: assertEquals("Nautical mile at North pole", 1861.67, e
102: .orthodromicDistance(0, 90 - 2 * hm, 0, 90), 0.2);
103: assertEquals("Nautical mile at South pole", 1861.67, e
104: .orthodromicDistance(0, 2 * hm - 90, 0, -90), 0.2);
105: assertEquals("International nautical mile", 1852.00, e
106: .orthodromicDistance(0, 45 - hm, 0, 45 + hm), 0.2);
107: for (double i = 0.01; i < 180; i += 1) {
108: final double base = 180 * random.nextDouble() - 90;
109: assertEquals(i + "° rotation", e.getSemiMajorAxis()
110: * Math.toRadians(i), e.orthodromicDistance(base, 0,
111: base + i, 0), 0.2);
112: }
113: /*
114: * Test the spherical model. The factory method should create
115: * a specialized class, which is not the usual Ellipsoid class.
116: */
117: final double radius = e.getSemiMajorAxis();
118: final double circumference = (radius * 1.00000001)
119: * (2 * Math.PI);
120: final DefaultEllipsoid s = DefaultEllipsoid.createEllipsoid(
121: "Sphere", radius, radius, e.getAxisUnit());
122: assertTrue("Spheroid class", !DefaultEllipsoid.class.equals(s
123: .getClass()));
124: for (double i = 0; i <= 180; i += 1) {
125: final double base = 360 * random.nextDouble() - 180;
126: assertEquals(i + "° rotation", s.getSemiMajorAxis()
127: * Math.toRadians(i), s.orthodromicDistance(base, 0,
128: base + i, 0), 0.001);
129: }
130: for (double i = -90; i <= +90; i += 1) {
131: final double meridian = 360 * random.nextDouble() - 180;
132: assertEquals(i + "° rotation", s.getSemiMajorAxis()
133: * Math.toRadians(Math.abs(i)), s
134: .orthodromicDistance(meridian, 0, meridian, i),
135: 0.001);
136: }
137: for (int i = 0; i < 100; i++) {
138: final double y1 = -90 + 180 * random.nextDouble();
139: final double y2 = -90 + 180 * random.nextDouble();
140: final double x1 = -180 + 360 * random.nextDouble();
141: final double x2 = -180 + 360 * random.nextDouble();
142: final double distance = s.orthodromicDistance(x1, y1, x2,
143: y2);
144: assertTrue("Range of legal values", distance >= 0
145: && distance <= circumference);
146: }
147: }
148:
149: /**
150: * Tests the {@link GeocentricTransform} class.
151: */
152: public void testGeocentricTransform() throws FactoryException,
153: TransformException {
154: /*
155: * Gets the math transform from WGS84 to a geocentric transform.
156: */
157: final DefaultEllipsoid ellipsoid = DefaultEllipsoid.WGS84;
158: final CoordinateReferenceSystem sourceCRS = DefaultGeographicCRS.WGS84_3D;
159: final CoordinateReferenceSystem targetCRS = DefaultGeocentricCRS.CARTESIAN;
160: final CoordinateOperation operation = opFactory
161: .createOperation(sourceCRS, targetCRS);
162: final MathTransform transform = operation.getMathTransform();
163: final int dimension = transform.getSourceDimensions();
164: assertEquals("Source dimension", 3, dimension);
165: assertEquals("Target dimension", 3, transform
166: .getTargetDimensions());
167: assertSame("Inverse transform", transform, transform.inverse()
168: .inverse());
169: assertInterfaced(transform);
170: /*
171: * Construct an array of 850 random points. The first 8 points
172: * are initialized to know values. Other points are left random.
173: */
174: final double cartesianDistance[] = new double[4];
175: final double orthodromicDistance[] = new double[4];
176: final double[] array0 = new double[900]; // Must be divisible by 3.
177: for (int i = 0; i < array0.length; i++) {
178: final int range;
179: switch (i % 3) {
180: case 0:
181: range = 360;
182: break; // Longitude
183: case 1:
184: range = 180;
185: break; // Latitidue
186: case 2:
187: range = 10000;
188: break; // Altitude
189: default:
190: range = 0;
191: break; // Should not happen
192: }
193: array0[i] = range * random.nextDouble() - (range / 2);
194: }
195: array0[0] = 35.0;
196: array0[1] = 24.0;
197: array0[2] = 8000; // 24°N 35°E 8km
198: array0[3] = 34.8;
199: array0[4] = 24.7;
200: array0[5] = 5000; // ... about 80 km away
201: cartesianDistance[0] = 80284.00;
202: orthodromicDistance[0] = 80302.99; // Not really exact.
203:
204: array0[6] = 0;
205: array0[7] = 0.0;
206: array0[8] = 0;
207: array0[9] = 180;
208: array0[10] = 0.0;
209: array0[11] = 0; // Antipodes; distance should be 2*6378.137 km
210: cartesianDistance[1] = ellipsoid.getSemiMajorAxis() * 2;
211: orthodromicDistance[1] = ellipsoid.getSemiMajorAxis() * Math.PI;
212:
213: array0[12] = 0;
214: array0[13] = -90;
215: array0[14] = 0;
216: array0[15] = 180;
217: array0[16] = +90;
218: array0[17] = 0; // Antipodes; distance should be 2*6356.752 km
219: cartesianDistance[2] = ellipsoid.getSemiMinorAxis() * 2;
220: orthodromicDistance[2] = 20003931.46;
221:
222: array0[18] = 95;
223: array0[19] = -38;
224: array0[20] = 0;
225: array0[21] = -85;
226: array0[22] = +38;
227: array0[23] = 0; // Antipodes
228: cartesianDistance[3] = 12740147.19;
229: orthodromicDistance[3] = 20003867.86;
230: /*
231: * Transform all points, and then inverse transform then. The resulting
232: * <code>array2</code> array should be equals to <code>array0</code>
233: * except for rounding errors. We tolerate maximal error of 0.1 second
234: * in longitude or latitude and 1 cm in height.
235: */
236: final double[] array1 = new double[array0.length];
237: final double[] array2 = new double[array0.length];
238: transform.transform(array0, 0, array1, 0, array0.length
239: / dimension);
240: transform.inverse().transform(array1, 0, array2, 0,
241: array1.length / dimension);
242: assertPointsEqual(
243: "transform(Geographic --> Geocentric --> Geographic)",
244: array0, array2, new double[] { 0.1 / 3600, 0.1 / 3600,
245: 0.01 });
246: /*
247: * Compare the distances between "special" points with expected distances.
248: * This test the ellipsoid orthodromic distance computation as well.
249: * We require a precision of 10 centimeters.
250: */
251: for (int i = 0; i < array0.length / 6; i++) {
252: final int base = i * 6;
253: final Point3d pt1 = new Point3d(array1[base + 0],
254: array1[base + 1], array1[base + 2]);
255: final Point3d pt2 = new Point3d(array1[base + 3],
256: array1[base + 4], array1[base + 5]);
257: final double cartesian = pt1.distance(pt2);
258: if (i < cartesianDistance.length) {
259: assertEquals("Cartesian distance[" + i + ']',
260: cartesianDistance[i], cartesian, 0.1);
261: }
262: /*
263: * Compare with orthodromic distance. Distance is computed using an ellipsoid
264: * at the maximal altitude (i.e. the length of semi-major axis is increased to
265: * fit the maximal altitude).
266: */
267: try {
268: final double altitude = Math.max(array0[base + 2],
269: array0[base + 5]);
270: final DefaultEllipsoid ellip = DefaultEllipsoid
271: .createFlattenedSphere("Temporary", ellipsoid
272: .getSemiMajorAxis()
273: + altitude, ellipsoid
274: .getInverseFlattening(), ellipsoid
275: .getAxisUnit());
276: double orthodromic = ellip.orthodromicDistance(
277: array0[base + 0], array0[base + 1],
278: array0[base + 3], array0[base + 4]);
279: orthodromic = XMath.hypot(orthodromic, array0[base + 2]
280: - array0[base + 5]);
281: if (i < orthodromicDistance.length) {
282: assertEquals("Orthodromic distance[" + i + ']',
283: orthodromicDistance[i], orthodromic, 0.1);
284: }
285: assertTrue("Distance consistency[" + i + ']',
286: cartesian <= orthodromic);
287: } catch (ArithmeticException exception) {
288: // Orthodromic distance computation didn't converge. Ignore...
289: }
290: }
291: }
292: }
|