sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Test TransverseMercator projection against data published in: Karney, C. F. F. (2009). Test data for the transverse Mercator projection [Data set]. Zenodo. Adjust the domain of validity threshold according the result of this test.
Date Sat, 25 May 2019 14:37:24 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new 25cec0a  Test TransverseMercator projection against data published in: Karney, C.
F. F. (2009). Test data for the transverse Mercator projection [Data set]. Zenodo. Adjust
the domain of validity threshold according the result of this test.
25cec0a is described below

commit 25cec0a6f08e0fdffc7ef723ae61ab161001c94b
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Sat May 25 16:35:48 2019 +0200

    Test TransverseMercator projection against data published in:
    Karney, C. F. F. (2009). Test data for the transverse Mercator projection [Data set].
Zenodo.
    Adjust the domain of validity threshold according the result of this test.
---
 .../operation/projection/TransverseMercator.java   | 43 +++++++++++++-
 .../projection/MapProjectionTestCase.java          |  1 +
 .../projection/TransverseMercatorTest.java         | 66 ++++++++++++++++++++--
 .../java/org/apache/sis/test/OptionalTestData.java | 27 +++++++++
 4 files changed, 130 insertions(+), 7 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
index a3067b3..cb6e414 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
@@ -74,6 +74,25 @@ public class TransverseMercator extends NormalizedProjection {
     private static final long serialVersionUID = -627685138188387835L;
 
     /**
+     * Distance from central meridian, in degrees, at which errors are considered too important.
+     * This threshold is determined by comparisons of computed values against values provided
by
+     * <cite>Karney (2009) Test data for the transverse Mercator projection</cite>
data file.
+     * On the WGS84 ellipsoid we observed:
+     *
+     * <ul>
+     *   <li>For ∆λ below 60°, errors below 1 centimetre.</li>
+     *   <li>For ∆λ between 60° and 66°, errors up to 0.1 metre.</li>
+     *   <li>For ∆λ between 66° and 70°, errors up to 1 metre.</li>
+     *   <li>For ∆λ between 70° and 72°, errors up to 2 metres.</li>
+     *   <li>For ∆λ between 72° and 74°, errors up to 10 metres.</li>
+     *   <li>For ∆λ between 74° and 76°, errors up to 30 metres.</li>
+     *   <li>For ∆λ between 76° and 78°, errors up to 1 kilometre.</li>
+     *   <li>For ∆λ greater than 85°, results are chaotic.</li>
+     * </ul>
+     */
+    static final double DOMAIN_OF_VALIDITY = 70;
+
+    /**
      * {@code false} for using the original formulas as published by EPSG, or {@code true}
for using formulas
      * modified using trigonometric identities. The use of trigonometric identities is for
reducing the amount
      * of calls to the {@link Math#sin(double)} and similar methods. Some identities used
are:
@@ -350,7 +369,7 @@ public class TransverseMercator extends NormalizedProjection {
                             final boolean derivate) throws ProjectionException
     {
         final double λ = srcPts[srcOff];
-        if (abs(λ) >= (0.9*PI/2)) {
+        if (abs(λ) >= DOMAIN_OF_VALIDITY * (PI/180)) {
             /*
              * The Transverse Mercator projection is conceptually a Mercator projection rotated
by 90°.
              * In Mercator projection, y values tend toward infinity for latitudes close
to ±90°.
@@ -652,9 +671,27 @@ public class TransverseMercator extends NormalizedProjection {
                                 final double[] dstPts, final int dstOff,
                                 final boolean derivate) throws ProjectionException
         {
-            final double λ = srcPts[srcOff  ];
+            final double λ = srcPts[srcOff];
+            /*
+             * The Transverse Mercator projection is conceptually a Mercator projection rotated
by 90°.
+             * In Mercator projection, y values tend toward infinity for latitudes close
to ±90° while
+             * in Transverse Mercator, x values tend toward infinity for longitudes close
to ±90° from
+             * central meridian (and at equator). After we pass the 90° limit, the Transverse
Mercator
+             * results at (90° + Δ) are the same as for (90° - Δ).
+             *
+             * Problem is that 90° is an ordinary longitude value, not even close to the
limit of longitude
+             * values range (±180°).  So having f(π/2+Δ, φ) = f(π/2-Δ, φ) results
in wrong behavior in some
+             * algorithms like the one used by Envelopes.transform(CoordinateOperation, Envelope).
 Since a
+             * distance of 90° from central meridian is outside projection domain of validity
anyway, we do
+             * not let the user go further.
+             *
+             * Note that in the case of ellipsoidal formulas (implemented by parent class),
we put the
+             * limit to a lower value (DOMAIN_OF_VALIDITY) because experience shows that
results close
+             * to equator become chaotic after 85° when using WGS84 ellipsoid. We do not
need to reduce
+             * the limit for the spherical formulas, because the mathematic are simpler and
the function
+             * still smooth until 90°.
+             */
             if (abs(λ) > PI/2) {
-                // See comment in the overridden class.
                 throw new ProjectionException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
             }
             final double φ    = srcPts[srcOff+1];
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
index 614c3b9..fc8d8c5 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
@@ -104,6 +104,7 @@ abstract strictfp class MapProjectionTestCase extends MathTransformTestCase
{
 
     /**
      * Initializes a complete projection (including conversion from degrees to radians) for
the given provider.
+     * Base CRS axis order is (longitude, latitude).
      */
     final void createCompleteProjection(final DefaultOperationMethod provider,
             final double semiMajor,
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
index 655058f..bc79bb3 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
@@ -16,27 +16,35 @@
  */
 package org.apache.sis.referencing.operation.projection;
 
+import java.util.Arrays;
 import java.util.Random;
+import java.io.IOException;
+import java.io.LineNumberReader;
 import org.opengis.util.FactoryException;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.internal.referencing.provider.TransverseMercatorSouth;
-import org.apache.sis.parameter.Parameters;
 import org.apache.sis.referencing.operation.transform.CoordinateDomain;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.util.CharSequences;
+import org.apache.sis.test.OptionalTestData;
 import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.DependsOn;
 import org.junit.Test;
 
 import static java.lang.Double.NaN;
+import static java.lang.StrictMath.min;
 import static java.lang.StrictMath.toRadians;
 import static org.apache.sis.test.Assert.*;
+import static org.apache.sis.internal.referencing.provider.TransverseMercator.LATITUDE_OF_ORIGIN;
+import org.opengis.test.CalculationType;
 
 
 /**
  * Tests the {@link TransverseMercator} class.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.7
+ * @version 1.0
  * @since   0.6
  * @module
  */
@@ -51,7 +59,7 @@ public final strictfp class TransverseMercatorTest extends MapProjectionTestCase
         final org.apache.sis.internal.referencing.provider.TransverseMercator method =
                 new org.apache.sis.internal.referencing.provider.TransverseMercator();
         final Parameters parameters = parameters(method, ellipsoidal);
-        parameters.getOrCreate(org.apache.sis.internal.referencing.provider.TransverseMercator.LATITUDE_OF_ORIGIN).setValue(latitudeOfOrigin);
+        parameters.getOrCreate(LATITUDE_OF_ORIGIN).setValue(latitudeOfOrigin);
         transform = new TransverseMercator(method, parameters);
         tolerance = NORMALIZED_TOLERANCE;
         validate();
@@ -103,7 +111,7 @@ public final strictfp class TransverseMercatorTest extends MapProjectionTestCase
                 0.5,        // Central meridian
                 2.5,        // Latitude of origin
                 NaN,        // Standard parallel 1 (none)
-                NaN,        // Standard parallel 1 (none)
+                NaN,        // Standard parallel 2 (none)
                 0.997,      // Scale factor
                 200,        // False easting
                 100);       // False northing
@@ -171,4 +179,54 @@ public final strictfp class TransverseMercatorTest extends MapProjectionTestCase
         tolerance = Formulas.LINEAR_TOLERANCE;
         verifyTransform(source, target);
     }
+
+    /**
+     * Compares with <cite>Karney (2009) Test data for the transverse Mercator projection</cite>.
+     * This is an optional test executed only if the {@code $SIS_DATA/Tests/TMcoords.dat}
file is found.
+     *
+     * @throws IOException if an error occurred while reading the test file.
+     * @throws FactoryException if an error occurred while creating the map projection.
+     * @throws TransformException if an error occurred while transforming coordinates.
+     */
+    @Test
+    public void compareAgainstDataset() throws IOException, FactoryException, TransformException
{
+        try (LineNumberReader reader = OptionalTestData.TRANSVERSE_MERCATOR.reader()) {
+            createCompleteProjection(new org.apache.sis.internal.referencing.provider.TransverseMercator(),
+                    WGS84_A,    // Semi-major axis length
+                    WGS84_B,    // Semi-minor axis length
+                    0,          // Central meridian
+                    0,          // Latitude of origin
+                    NaN,        // Standard parallel 1 (none)
+                    NaN,        // Standard parallel 2 (none)
+                    0.9996,     // Scale factor
+                    0,          // False easting
+                    0);         // False northing
+            final double[] source = new double[2];
+            final double[] target = new double[2];
+            String line;
+            while ((line = reader.readLine()) != null) {
+                Arrays.fill(source, Double.NaN);
+                final CharSequence[] split = CharSequences.split(line, ' ');
+                for (int i=min(split.length, 4); --i >= 0;) {
+                    final double value = Double.parseDouble(split[i].toString());
+                    if (i <= 1) source[i ^ 1] = value;                              //
Swap axis order.
+                    else        target[i - 2] = value;
+                }
+                // Relax tolerance for longitudes very far from central meridian.
+                final double longitude = StrictMath.abs(source[0]);
+                if (longitude < TransverseMercator.DOMAIN_OF_VALIDITY) {
+                    if (StrictMath.abs(source[1]) >= 89.9) tolerance = 0.1;
+                    else if (longitude <= 60) tolerance = Formulas.LINEAR_TOLERANCE;
+                    else if (longitude <= 66) tolerance = 0.1;
+                    else if (longitude <= 70) tolerance = 1;
+                    else if (longitude <= 72) tolerance = 2;
+                    else if (longitude <= 74) tolerance = 10;
+                    else if (longitude <= 76) tolerance = 30;
+                    else                      tolerance = 1000;
+                    transform.transform(source, 0, source, 0, 1);
+                    assertCoordinateEquals(line, target, source, reader.getLineNumber(),
CalculationType.DIRECT_TRANSFORM);
+                }
+            }
+        }
+    }
 }
diff --git a/core/sis-utility/src/test/java/org/apache/sis/test/OptionalTestData.java b/core/sis-utility/src/test/java/org/apache/sis/test/OptionalTestData.java
index 5abdad4..176a297 100644
--- a/core/sis-utility/src/test/java/org/apache/sis/test/OptionalTestData.java
+++ b/core/sis-utility/src/test/java/org/apache/sis/test/OptionalTestData.java
@@ -39,6 +39,33 @@ import static org.junit.Assume.assumeNotNull;
  */
 public enum OptionalTestData {
     /**
+     * Transverse Mercator projection on the WGS84 ellipsoid. Central meridian is 0° and
scale factor is 0.9996.
+     *
+     * <dl>
+     *   <dt>File:</dt>
+     *   <dd>{@code TMcoords.dat}</dd>
+     *   <dt>Size:</dt>
+     *   <dd>34422206 bytes (34 Mb)</dd>
+     *   <dt>MD5 sum:</dt>
+     *   <dd>{@code 91b817eae34aef3cd8d67b5a4e8e798f}</dd>
+     *   <dt>Source:</dt>
+     *   <dd><a href="http://doi.org/10.5281/zenodo.32470">Karney, C. F. F. (2009).
Test data for the transverse Mercator projection [Data set]. Zenodo.</a></dd>
+     * </dl>
+     *
+     * Each line in the test file gives the following numbers (space delimited):
+     *
+     * <ol>
+     *   <li>φ — latitude (degrees)</li>
+     *   <li>λ — longitude (degrees)</li>
+     *   <li>E — easting (metres)</li>
+     *   <li>N — Northing (metres)</li>
+     *   <li>meridian convergence (degrees)</li>
+     *   <li>scale</li>
+     * </ol>
+     */
+    TRANSVERSE_MERCATOR("TMcoords.dat"),
+
+    /**
      * Geodesic distances, rhumb line length and azimuths on WGS84 ellipsoid computed from
a set of points.
      *
      * <dl>


Mime
View raw message