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
commit 440c700f6b58223b3bb5c074da3ddcf2e3ee2c85
Author: Martin Desruisseaux <martin.desruisse...@geomatys.com>
AuthorDate: Sat Jan 8 07:22:04 2022 +0100
Fix a NaN value caused by rounding error in Oblique Mercator inverse
projection.
https://issues.apache.org/jira/browse/SIS-532
---
.../operation/projection/ObliqueMercator.java | 11 ++++++--
.../operation/projection/ObliqueMercatorTest.java | 31 +++++++++++++++++++++-
2 files changed, 39 insertions(+), 3 deletions(-)
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java
index c4380df..7203e1f 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java
@@ -55,7 +55,7 @@ import static
org.apache.sis.internal.referencing.provider.ObliqueMercatorCenter
* @author Rueben Schulz (UBC)
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.0
+ * @version 1.2
*
* @see Mercator
* @see TransverseMercator
@@ -393,6 +393,13 @@ public class ObliqueMercator extends ConformalProjection {
* + sin(8χ)⋅(e⁸⋅4279/161280)
*
* The calculation of χ and its use in serie expansion is performed by
the φ(t) function.
+ *
+ * Note: U shall be in the [−1 … +1] range, otherwise we get φ = NaN.
At the limit U = ±1 we get φ = ±π/2
+ * (division by zero is okay: infinite intermediate values produce the
correct final result with φ(∞)=-π/2).
+ * Empirical tests suggest that U is inside valid range even with
points far outside the projection domain.
+ * We found out-of-range values only at short distances (about 0.01
meter) from pole. In those cases U was
+ * out-of-range only by 1 ULP. The (abs(U) >= 1) check is for handling
those cases. Note that if U is NaN,
+ * then we want to use the general formula for propagating the NaN.
*/
final double x = srcPts[srcOff ];
final double y = srcPts[srcOff+1];
@@ -403,7 +410,7 @@ public class ObliqueMercator extends ConformalProjection {
final double V = sin(y);
final double U = (V*cosγ0 + S*sinγ0) / T;
final double λ = -atan2((S*cosγ0 - V*sinγ0), cos(y)) / B;
- final double φ = φ(pow(H / sqrt((1 + U) / (1 - U)), 1/B));
+ final double φ = (abs(U) >= 1) ? copySign(PI/2, U) : φ(pow(H /
sqrt((1 + U) / (1 - U)), 1/B));
dstPts[dstOff ] = λ;
dstPts[dstOff+1] = φ;
}
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueMercatorTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueMercatorTest.java
index 67e8d70..a3318ea 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueMercatorTest.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueMercatorTest.java
@@ -16,11 +16,14 @@
*/
package org.apache.sis.referencing.operation.projection;
+import org.opengis.test.ToleranceModifier;
import org.opengis.util.FactoryException;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.referencing.datum.Ellipsoid;
import org.opengis.referencing.operation.TransformException;
+import org.opengis.referencing.operation.MathTransformFactory;
import org.apache.sis.internal.referencing.provider.ObliqueMercatorCenter;
+import org.apache.sis.internal.system.DefaultFactories;
import org.apache.sis.referencing.CommonCRS;
import org.apache.sis.parameter.Parameters;
import org.apache.sis.test.DependsOn;
@@ -34,7 +37,7 @@ import static java.lang.StrictMath.*;
*
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.0
+ * @version 1.2
* @since 1.0
* @module
*/
@@ -96,6 +99,32 @@ public final strictfp class ObliqueMercatorTest extends
MapProjectionTestCase {
}
/**
+ * Tests with a latitude close to 90°.
+ *
+ * @throws FactoryException if an error occurred while creating the map
projection.
+ * @throws TransformException if an error occurred while converting a
point.
+ *
+ * <a href="https://issues.apache.org/jira/browse/SIS-532";>SIS-532</a>
+ */
+ @Test
+ public void testPole() throws TransformException, FactoryException {
+ tolerance = 0.01;
+ transform = create(179.8, 89.8,
-174).createMapProjection(DefaultFactories.forBuildin(MathTransformFactory.class));
+ transform = transform.inverse();
+ validate();
+ /*
+ * The projection of (180, 90) with SIS 1.1 is (+0.004715980030596256,
22338.795490272343).
+ * Empirical cordinates shifted by 0.01 meter: (-0.005463426921067797,
22338.792057282844).
+ * With those shifted coordinated, Apache SIS 1.1 was used to compute
φ = NaN because the
+ * U′ value in `ObliqueMercator.inverseTransform(…)` was slightly
greater than 1.
+ */
+ isInverseTransformSupported = false;
+ toleranceModifier = ToleranceModifier.GEOGRAPHIC;
+ verifyTransform(new double[] {-0.005464, 22338.792057},
+ new double[] {300, 90});
+ }
+
+ /**
* Tests the <cite>"Hotine Oblique Mercator (variant B)"</cite>
(EPSG:9815) projection method.
* This test is defined in GeoAPI conformance test suite.
*
