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: Finer determination of Transverse Mercator domain of validity.
Date Wed, 24 Feb 2021 15:13:57 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 27a08d6  Finer determination of Transverse Mercator domain of validity.
27a08d6 is described below

commit 27a08d6d693d3bbd87b7e92f8cb714da6f5bd258
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Feb 24 16:12:51 2021 +0100

    Finer determination of Transverse Mercator domain of validity.
---
 .../operation/projection/TransverseMercator.java   |  39 +++++++++++++++------
 .../projection/TransverseMercatorTest.java         |  22 ++++++++----
 .../doc-files/TransverseMercatorErrors.png         | Bin 0 -> 7954 bytes
 3 files changed, 44 insertions(+), 17 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 837c488..fb00694 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
@@ -59,7 +59,7 @@ import static org.apache.sis.internal.referencing.provider.TransverseMercator.*;
  *
  * @author  Martin Desruisseaux (Geomatys)
  * @author  Rémi Maréchal (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @see Mercator
  * @see ObliqueMercator
@@ -77,7 +77,7 @@ public class TransverseMercator extends NormalizedProjection {
      * 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:
+     * On the WGS84 ellipsoid we observed at equator:
      *
      * <ul>
      *   <li>For ∆λ below 60°, errors below 1 centimetre.</li>
@@ -89,8 +89,24 @@ public class TransverseMercator extends NormalizedProjection {
      *   <li>For ∆λ between 76° and 78°, errors up to 1 kilometre.</li>
      *   <li>For ∆λ greater than 85°, results are chaotic.</li>
      * </ul>
+     *
+     * For latitudes greater than 20°, errors are less than 0.5 meters for all ∆λ &lt;
83°.
+     * Errors become suddenly much higher at ∆λ &gt; 82.6° no matter the latitude.
+     */
+    static final double DOMAIN_OF_VALIDITY = 82.5;
+
+    /**
+     * The domain of validity at equator. We keep using the limit at all latitudes up to
+     * {@value #LATITUDE_OF_REDUCED_DOMAIN} degrees, even if actually the limit could be
+     * progressively relaxed until it reaches {@value #DOMAIN_OF_VALIDITY} degrees.
      */
-    static final double DOMAIN_OF_VALIDITY = 70;
+    static final double DOMAIN_OF_VALIDITY_AT_EQUATOR = 70;
+
+    /**
+     * The maximal latitude (exclusive) where to use {@link #DOMAIN_OF_VALIDITY_AT_EQUATOR}
+     * instead of {@link #DOMAIN_OF_VALIDITY}.
+     */
+    static final double LATITUDE_OF_REDUCED_DOMAIN = 20;
 
     /**
      * {@code false} for using the original formulas as published by EPSG, or {@code true}
for using formulas
@@ -376,7 +392,8 @@ public class TransverseMercator extends NormalizedProjection {
                             final boolean derivate) throws ProjectionException
     {
         final double λ = srcPts[srcOff];
-        if (abs(λ) >= DOMAIN_OF_VALIDITY * (PI/180)) {
+        final double φ = srcPts[srcOff+1];
+        if (abs(λ) >= DOMAIN_OF_VALIDITY_AT_EQUATOR * (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°.
@@ -395,11 +412,13 @@ public class TransverseMercator extends NormalizedProjection {
              * 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 (Math.abs(IEEEremainder(λ, 2*PI)) >= DOMAIN_OF_VALIDITY * (PI/180)) {
   // More costly check.
+            final double limit = (abs(φ) < LATITUDE_OF_REDUCED_DOMAIN * (PI/180))
+                                    ? (PI/180) * DOMAIN_OF_VALIDITY_AT_EQUATOR
+                                    : (PI/180) * DOMAIN_OF_VALIDITY;
+            if (Math.abs(IEEEremainder(λ, 2*PI)) >= limit) {
                 throw new ProjectionException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
             }
         }
-        final double φ     = srcPts[srcOff+1];
         final double sinλ  = sin(λ);
         final double ℯsinφ = sin(φ) * eccentricity;
         final double Q     = asinh(tan(φ)) - atanh(ℯsinφ) * eccentricity;
@@ -681,6 +700,9 @@ public class TransverseMercator extends NormalizedProjection {
                                 final boolean derivate) throws ProjectionException
         {
             final double λ = srcPts[srcOff];
+            final double φ = srcPts[srcOff + 1];
+            final double sinλ = sin(λ);
+            final double cosλ = cos(λ);
             /*
              * 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
@@ -700,12 +722,9 @@ public class TransverseMercator extends NormalizedProjection {
              * the limit for the spherical formulas, because the mathematic are simpler and
the function
              * still smooth until 90°.
              */
-            if (abs(λ) > PI/2) {
+            if (cosλ < 0) {                     // Implies Math.abs(IEEEremainder(λ,
2*PI)) > PI/2
                 throw new ProjectionException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
             }
-            final double φ    = srcPts[srcOff+1];
-            final double sinλ = sin(λ);
-            final double cosλ = cos(λ);
             final double sinφ = sin(φ);
             final double cosφ = cos(φ);
             final double tanφ = sinφ / cosφ;
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 164d1ac..e526e3f 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,7 +16,6 @@
  */
 package org.apache.sis.referencing.operation.projection;
 
-import java.util.Arrays;
 import java.util.Random;
 import java.io.IOException;
 import java.io.LineNumberReader;
@@ -33,7 +32,7 @@ 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.abs;
 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;
@@ -183,6 +182,12 @@ public final strictfp class TransverseMercatorTest extends MapProjectionTestCase
     /**
      * 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.
+     * The errors observed as of February 2021 are illustrated below. In this image of size
91×91 pixels,
+     * pixels coordinates are longitude and latitude coordinates with (0,0) in the upper-left
corner.
+     * Black pixels are errors less than 0.5 meters. Yellow pixels are errors between 0.5
and 1.5 meters.
+     * Red pixels are errors of 254.5 meters or more.
+     *
+     * <img src="doc-files/TransverseMercatorErrors.png" alt="Errors of Transverse Mercator
projection">
      *
      * @throws IOException if an error occurred while reading the test file.
      * @throws FactoryException if an error occurred while creating the map projection.
@@ -205,17 +210,20 @@ public final strictfp class TransverseMercatorTest extends MapProjectionTestCase
             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;) {
+                for (int i=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;
+                final double longitude = abs(source[0]);
+                final double latitude  = abs(source[1]);
+                final double limit     = (latitude < TransverseMercator.LATITUDE_OF_REDUCED_DOMAIN
+                                                   ? TransverseMercator.DOMAIN_OF_VALIDITY_AT_EQUATOR
+                                                   : TransverseMercator.DOMAIN_OF_VALIDITY);
+                if (longitude < limit) {
+                    if (latitude >= 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;
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/doc-files/TransverseMercatorErrors.png
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/doc-files/TransverseMercatorErrors.png
new file mode 100644
index 0000000..47a76b7
Binary files /dev/null and b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/doc-files/TransverseMercatorErrors.png
differ


Mime
View raw message