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: Implement ellipsoidal formulas for the sinusoidale projection. This commit completes https://issues.apache.org/jira/browse/SIS-450
Date Thu, 25 Apr 2019 13:39:25 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 becfcc7  Implement ellipsoidal formulas for the sinusoidale projection. This commit
completes https://issues.apache.org/jira/browse/SIS-450
becfcc7 is described below

commit becfcc7847aad755d1cb1156eeac7b76b0625902
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Thu Apr 25 15:38:27 2019 +0200

    Implement ellipsoidal formulas for the sinusoidale projection.
    This commit completes https://issues.apache.org/jira/browse/SIS-450
---
 .../operation/projection/MeridianArcBased.java     |  3 +-
 .../operation/projection/Sinusoidal.java           | 43 +++++++++++++++------
 .../operation/projection/SinusoidalTest.java       | 45 ++++++++++++++++------
 3 files changed, 65 insertions(+), 26 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
index c624aa0..db20419 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
@@ -230,9 +230,8 @@ abstract class MeridianArcBased extends NormalizedProjection {
     final double latitude(final double distance) {
         double    φ  = distance / rld;            // Rectifying latitude µ used as first
approximation.
         double sinφ  = sin(φ);
-        double cosφ  = cos(φ);
         double sinφ2 = sinφ * sinφ;
-        φ += sinφ*cosφ*(ci1 + sinφ2*(ci2 + sinφ2*(ci3 + sinφ2*ci4)));             
     // Snyder 3-26.
+        φ += cos(φ)*sinφ*(ci1 + sinφ2*(ci2 + sinφ2*(ci3 + sinφ2*ci4)));           
     // Snyder 3-26.
         /*
          * We could improve accuracy by continuing from here with Newton's iterative method
          * (see MeridianArcTest.inverse(…) for implementation). However those iterations
requires
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Sinusoidal.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Sinusoidal.java
index dd22af5..6180442 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Sinusoidal.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Sinusoidal.java
@@ -118,15 +118,28 @@ public class Sinusoidal extends MeridianArcBased {
                             final double[] dstPts, final int dstOff,
                             final boolean derivate) throws ProjectionException
     {
-        final double λ    = srcPts[srcOff  ];
-        final double φ    = srcPts[srcOff+1];
-        final double cosφ = cos(φ);
-        // TODO: replace by ellipsoidal formulas.
+        final double λ     = srcPts[srcOff  ];
+        final double φ     = srcPts[srcOff+1];
+        final double cosφ  = cos(φ);
+        final double sinφ  = sin(φ);
+        final double sinφ2 = sinφ * sinφ;
+        final double rν2   = 1 - eccentricitySquared * sinφ2;
+        final double rν    = sqrt(rν2);                                 // Reciprocal of
the radius of curvature.
+        final double dx_dλ = cosφ / rν;                                 // Part of Snyder
30-8.
+        /*
+         * Note: in theory x/cos(φ) is indeterminate at φ=±π/2. However in this code,
+         * that indetermination never happen because there is no exact representation
+         * of π/2 in base 2, so cos(φ) can never return 0.
+         */
         if (dstPts != null) {
-            dstPts[dstOff]   = λ * cosφ;
-            dstPts[dstOff+1] = φ;
+            dstPts[dstOff  ] = dx_dλ * λ;
+            dstPts[dstOff+1] = distance(φ, sinφ, cosφ);
+        }
+        if (!derivate) {
+            return null;
         }
-        return derivate ? new Matrix2(cosφ, -λ*sin(φ), 0, 1) : null;
+        final double dx_dφ = λ * sinφ * (eccentricitySquared*(cosφ*cosφ) / rν2 - 1)
/ rν;
+        return new Matrix2(dx_dλ, dx_dφ, 0, dM_dφ(sinφ2));
     }
 
     /**
@@ -137,11 +150,12 @@ public class Sinusoidal extends MeridianArcBased {
     protected void inverseTransform(final double[] srcPts, final int srcOff,
                                     final double[] dstPts, final int dstOff)
     {
-        final double x = srcPts[srcOff  ];
-        final double φ = srcPts[srcOff+1];
-        // TODO: replace by ellipsoidal formulas.
-        dstPts[dstOff  ] = x / cos(φ);
-        dstPts[dstOff+1] = φ;
+        final double x    = srcPts[srcOff  ];
+        final double y    = srcPts[srcOff+1];
+        final double φ    = latitude(y);
+        final double sinφ = sin(φ);
+        dstPts[dstOff  ]  = x * sqrt(1 - eccentricitySquared * (sinφ*sinφ)) / cos(φ);
          // Part of Snyder 30-11
+        dstPts[dstOff+1]  = φ;
     }
 
 
@@ -214,6 +228,11 @@ public class Sinusoidal extends MeridianArcBased {
         {
             final double x = srcPts[srcOff  ];
             final double φ = srcPts[srcOff+1];
+            /*
+             * Note: in theory x/cos(φ) is indeterminate at φ=±π/2. However in this code,
+             * that indetermination never happen because there is no exact representation
+             * of π/2 in base 2, so cos(φ) can never return 0.
+             */
             dstPts[dstOff  ] = x / cos(φ);              // Part of Snyder 30-5
             dstPts[dstOff+1] = φ;                       // Part of Snyder 30-6
         }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
index ae0f141..a448db2 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
@@ -61,11 +61,13 @@ public final strictfp class SinusoidalTest extends MapProjectionTestCase
{
         createProjection(false);
         assertTrue(isInverseTransformSupported);
         verifyTransform(
-            new double[] {              // (λ,φ) coordinates in degrees to project.
-                  2,      1
+            new double[] {                  // (λ,φ) coordinates in degrees to project.
+                  2,              1,
+                -75 - -90,      -50         // Snyder example is relative to λ₀ = 90°W.
             },
-            new double[] {              // Expected (x,y) results in metres.
-               223368.12, 111701.07     // Values taken from PROJ.4.
+            new double[] {                  // Expected (x,y) results in metres.
+               223368.12,    111701.07,     // Values taken from PROJ.4.
+              1077000.98,  -5585053.61      // Values derived from Snyder page 365.
             });
     }
 
@@ -76,29 +78,30 @@ public final strictfp class SinusoidalTest extends MapProjectionTestCase
{
      * @throws TransformException if an error occurred while projecting a point.
      */
     @Test
-    @org.junit.Ignore("Ellipsoidal formula not yet implemented.")
     public void testEllipsoidal() throws FactoryException, TransformException {
         createProjection(true);
         assertTrue(isInverseTransformSupported);
         verifyTransform(
-            new double[] {              // (λ,φ) coordinates in degrees to project.
-                  2,      1
+            new double[] {                  // (λ,φ) coordinates in degrees to project.
+                  2,              1,
+                -75 - -90,      -50         // Snyder example is relative to λ₀ = 90°W.
             },
-            new double[] {              // Expected (x,y) results in metres.
-               222605.30, 110574.39     // Values taken from PROJ.4.
+            new double[] {                  // Expected (x,y) results in metres.
+               222605.30,    110574.39,     // Values taken from PROJ.4.
+              1075436.30,  -5540847.04      // Snyder values modified for WGS84 (Δx ≈
35m and Δy ≈ 219m).
             });
     }
 
     /**
-     * Tests the derivatives at a few points. This method compares the derivatives computed
by
-     * the projection with an estimation of derivatives computed by the finite differences
method.
+     * Tests the derivatives at a few points on a sphere. This method compares the derivatives
computed
+     * by the projection with an estimation of derivatives computed by the finite differences
method.
      *
      * @throws FactoryException if an error occurred while creating the map projection.
      * @throws TransformException if an error occurred while projecting a point.
      */
     @Test
     @DependsOnMethod("testInverseDerivative")
-    public void testDerivative() throws FactoryException, TransformException {
+    public void testDerivativeOnSphere() throws FactoryException, TransformException {
         createProjection(false);
         final double delta = (100.0 / 60) / 1852;               // Approximatively 100 metres.
         derivativeDeltas = new double[] {delta, delta};
@@ -106,4 +109,22 @@ public final strictfp class SinusoidalTest extends MapProjectionTestCase
{
         verifyDerivative(15,  30);
         verifyDerivative(10, -60);
     }
+
+    /**
+     * Tests the derivatives at a few points on an ellipsoid. This method compares the derivatives
computed
+     * by the projection with an estimation of derivatives computed by the finite differences
method.
+     *
+     * @throws FactoryException if an error occurred while creating the map projection.
+     * @throws TransformException if an error occurred while projecting a point.
+     */
+    @Test
+    @DependsOnMethod("testInverseDerivative")
+    public void testDerivativeOnEllipsoid() throws FactoryException, TransformException {
+        createProjection(true);
+        final double delta = (100.0 / 60) / 1852;               // Approximatively 100 metres.
+        derivativeDeltas = new double[] {delta, delta};
+        tolerance = 1E-6;                                       // More severe than Formulas.LINEAR_TOLERANCE.
+        verifyDerivative(15,  30);
+        verifyDerivative(10, -60);
+    }
 }


Mime
View raw message